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CONTRACTOR REPORT 


VELOCITY-PRESSURE INTEGRATED VERSUS PENALTY FINITE ELEMENT 
METHODS FOR HIGH REYNOLDS NUMBER FLOWS 

I. INTRODUCTION 


Various finite element methods for the Navier-Stokes equations have been pro- 
posed during the last decades. These finite element methods may be categorized into 
three classes based on the way pressure has been treated; these are the velocity- 
pressure integrated mixed interpolation methods [1-3], the penalty methods [3-5], 
and the velocity-pressure segregated methods [6-8]. 

The velocity-pressure integrated, mixed interpolation methods do not require 
any approximation at the differential equation level. Whereas simplified pressure 
and/or pressure correction equations are used in the velocity-pressure segregated 
methods, and the penalized conservation of mass equation is used in the penalty methods. 
Conceptually, the velocity-pressure integrated methods would satisfy the conservation 
of mass equation most rigorously in a sense that no approximation has to be made at 
the level of differential equation. The most classical velocity-pressure integrated 
method is based on an 8-node velocity, 4-node pressure flow element. Unfortunately, 
this element yielded inaccurate pressure which became more severe as the Reynolds 
number was increased. It has been shown in Reference 2 that a 9-node velocity, 
3-node pressure flow element, when used in the velocity- pressure integrated finite 
element method, yielded accurate velocity and pressure for high Reynolds number 
flows. It has also been shown that no upwinding technique was necessary to obtain 
computational results which were free of numerical wiggles for high Reynolds number 
flows. 

The velocity -pres sure segregated methods have been motivated by the success 
of the finite difference computational methods based on segregated formulation of the 
Navier-Stokes equations, such as the SIMPLE (Semi- Implicit Method for Pressure- 
Linked Equations) algorithm [9]. The computational results obtained by using the 
segregated finite element methods have not shown, as yet, any significant improve- 
ment in accuracy compared with the other finite element methods. However, signifi- 
cant computational efficiency in computer time and storage can be achieved by the 
velocity- pressure segregated methods. 

In the penalty method, the pressure is pre- eliminated from the Navier-Stokes 
equations by penalizing the conservation of mass equation. The conservation of mass 
constraint could be satisfied rigorously as the penalty number approaches infinity. 
However, the frequently used penalty number has been limited to 10® through 10-1-0 
times of the kinematic viscosity in order to avoid ill-conditioned matrix [10]. The 
influence of the penalty number on the converged solution can be found in References 
3 to 5, among many others. The consistent penalty finite element method [4] has 
been used in the present study. The improvements realized by using the new pres- 
sure interpolation polynomials in the consistent penalty finite element method is dis- 
cussed in detail. It is shown that the consistent penalty method and the velocity- 
pressure integrated method yielded comparable computational -results in accuracy and 
convergence rate. 



The nonlinear, finite element system of equations has been solved by the direct 
(Picard) iteration method using a frontal solver [1,10]. It is intended to extend the 
present finite element code to solve turbulent flows. For turbulent flows, a strongly 
convergent solution technique, which may require severe under- relaxation, need to 
be used to obtain convergent solutions [11-13]. Thus, inclusion of the Newton- 
Raphson method into the present finite element code has not been considered. 


II. FINITE ELEMENT EQUATIONS 


A finite element system of equations for two-dimensional, laminar, steady, 
incompressible flows is described below. The method is based on the standard 
Gailerkin finite element method. In the following discussion, repeated indices imply 
summation over the indices, unless otherwise specified. 


The Navier-Stokes equations are given as: 



where ft is the open bounded domain, the subscripts i and j denote the coordinate 
directions, p is the density, Uj is the velocity component in the i-th coordinate 

direction, p is the pressure, y is the molecular viscosity of the fluid, b. is the body 

force in the i-th coordinate direction, and 6 .. is the Kronecker delta such that 

6^=1 for i = j and 6 . ^ = 0 for i ^ j . The boundary conditions are given as : 



for x €. 9 ft . 
~ 1 


T. = T..n. for x e 3ft 
l i] j 


(3) 




where x = (x,y), 3ft^ is part of the boundary on which Dirichlet boundary condition 
is specified, 3ftg is the rest of the boundary on which Newmann boundary condition 
is specified, T^ is the surface traction, and is the stress tensor given as ij. = 
y ( suj/axj+auj/axj ) - ps^. 

In the penalty method, the conservation of mass equation is expressed as: 
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where A is the penalty number. The finite element system of equations for the con- 
sistent penalty method is described below. 

The finite element system of equations has been obtained by the standard 
Galerkin method [14]. In the method: the flow domain is discretized into a number 

of elements; the Navier-Stokes equations are multiplied by appropriate test functions; 
the second order stress tensor term is integrated by parts using the Green- Gauss 
theorem; the continuous flow variables are interpolated using the nodal values of 
these variables and the appropriate interpolation polynomials; the weak form Navier- 
Stokes equations are integrated over each element to obtain the discrete, element 
system of equations; and the element system of equations are assembled to obtain the 
global system of equations. Detailed derivation of the finite element system of 
equations can be found in References 1, 2, and 15, among many others. 

The system of equations for an element (ft )'is given as, in matrix form: 
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u. is a column vector of nodal values of the velocity component u., p is a column 

^ -1 • *>_/ 

vector of nodal pressure, $ is a column vector of interpolating polynomials for velo- 
city, \p is a column vector of interpolating polynomials for pressure, {b.c.} is a 

t — ’ 

column vector of the flux boundary condition, and the subscripts i and j denote the 
spatial dimensions. The integrations in equations (7) thorugh (12) have been 
evaluated using the Gauss numerical quadrature method with three Gauss points in 
each coordinate direction. 

In the velocity-pressure integrated method, the right hand side of equation (6) 
has been replaced by a null column vector, and the element system of equations given 
as equations (5) and (6) were assembled to obtain the global system of equations. 

In the penalty finite element method, equation (6) has been inverted to obtain a 
column vector of the nodal pressure and the result has been substituted into equa- 
tion (5) to obtain: 
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The flow element used in the present study is briefly described below . The 
velocities were interpolated using the bi- quadratic shape functions and the pressure 
was interpolated using the linear shape functions defined on a triangular element. 
The triangular pressure element is contained inside the quadratic isoparametric ele- 
ment [ 23 . Tiie three pressure nodes are located at the three Gauss points of the 
three-point Gauss quadrature rule for -quadrilateral elements [16]. The coordinates 
of the pressure nodes on the computational element are given as: 
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where C n = (C n> n n ), and n denotes the pressure node number. The shape function 
for each pressure node is given as: 
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Slightly different pressure interpolation methods have also been tested in the 
present study. These pressure interpolation polynomials are given as [4]: 


/ = (1, x, y) (16) 

r-J 

and 

/ = (1, 5, n) . (17) 

These three sets of pressure interpolation polynomials, equations (15) to (17), 
belong to the same approximation space if rectangular elements are used; and equa- 
tions (15) and (17) belong to the same approximation space if arbitrary distorted 
quadrilateral elements are used. Thus, any pressure interpolation polynomials which 
belong to the same approximation space should yield identical computational results. 

V Any difference in the computational results has to be related to the matrix condition 

and the computer round off error [ 2] . The performance of these three sets of 
slightly different pressure interpolation polynomials, when used in the consistent 
penalty method, are discussed in the following section. 

The nonlinear system of equations has been solved by a direct (Picard) itera- 
tion method using a frontal solver. The solutions have been updated using an under- 
relaxation method given as: 


a.* = a a^ n + (1 - a) a. n_1 


(18) 
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where represents any degree of freedom, a is the under- relaxation number, the 
superscripts n and n-1 denote the iteration levels, and aj* is the updated solution, 
a = 0.8 and a = 1 have been used for the velocities and the pressure, respectively. 


III. EXAMPLE PROBLEMS 


The finite element methods described in the previous section were tested by 
solving a lid-driven cavity flow [15, 17-21], a laminar backward- facing step flow 
[22-24], a laminar flow through a nest of cylinders [25-27], and a channel flow with 
an internal blockage. For the cavity flow, sharp boundary layers develop along all 
the boundary edges of the cavity at high Reynolds numbers. For the backward- 
facing step flow, the flow expands abruptly at the convex corner of the step. These 
flows contain subtle pressure driven recirculation zones. Inside the pressure driven 
recirculation zones, the local Reynolds number may become vanishingly small. Due 
to these reasons, obtaining convergent solutions with any iterative numerical method 
can be quite difficult [23]. Therefore, these two flows provide serious test cases 
for any numerical methods. The laminar flow through a nest of cylinders has been 
considered to investigate the convergence nature of the two finite element methods. 
The channel flow with an internal blockage has been included herein to investigate 
the source of numerical wiggles for high Reynolds number flows. 

The error norm ( e^) has been defined as: 


Max 
e i “ j=l.,N 



(19) 
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where the subscript i (i = u, v, or p) denotes each component of the flow variables, 

a. . denotes the nodal value of the i-th flow variable for the j-th node, a. denotes 
i,j i,max 

the maximum value of the i-th flow variable in the previous iteration level, and N is 

the total number of nodes. Solving the coupled system of equations once has been 

counted as an iteration. 

In the penalty method, the pressure is recovered in the post process using the 
penalized conservation of mass equation. The quality of the recovered pressure 
depends on the velocity. In the present study, the pressure has been recovered at 
the end of each iteration. The purpose has been to provide some insight into the 
convergence nature of the pressure in the penalty methods. This information may be 
helpful in selecting the convergence criterion in application situations. Note that the 
required number of iterations to obtain a convergent solution depends on the pre- 
scribed convergence criterion in nonlinear problems. 

The pressure is discontinuous across element boundaries. Thus, the nodal 
pressure at the velocity node has been obtained by averaging all the pressure con- 
tributions made by the: elements containing the node; and each contribution has been 
evaluated using an equation given as: 


( 20 ) 
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where ^ is the pressure interpolation polynomial and pj is the corresponding pressure 
degree of freedom. 

A penalty number of (y/p) x 10^ has been used in the present study. 


3.1 Lid-Driven Cavity Flow 

A lid-driven cavity flow for Reynolds number of 10,000 is considered below. 

The fine grid finite difference computational results of the flow can be found in 
References 18 and 19. The no slip boundary condition (u = v = 0) has been applied 
at all the boundaries except at y = 1 where u = 1 and v = 0. A fixed pressure 
boundary condition has been prescribed at an arbitrary pressure node inside the 
domain. The Reynolds number is defined as R g = pUL/y , where U = 1 is the velocity 

of the lid, L = 1 is the reference length, and y is the molecular viscosity of the 
fluid. The computational domain has been discretized by unequally spaced 32 x 32 
quadratic elements [2]. The trivial solution (u = v = p = 0) has been used as an 
initial guess for all the cases. 

The streamline contours and the normalized pressure contours obtained by using 
the penalty method with the pressure interpolation polynomials given as equations ( 15) 
and (16) are shown in Figures 1 and 2, respectively. The normalized pressure (P) 
has been obtained from the static pressure (p) using a relationship given as [20] : 


P = pL/U/.y. 


( 21 ) 


The streamline and the pressure contour labels are given in Table 1. In Figures 1 
and 2, the reference pressure of p = 0 has been assigned at the middle of the bottom 
wall. 


TABLE 1. CONTOUR LABEL FOR CAVITY FLOW 


(a) Streamline Contour 


Label 


Label 

i> 

Label 


A 

-0.11 

F 

-0.03 

K 

2.X10' 4 

B 

-0.10 

G 

-0.01 

L 

5.X10' 4 

C 

-0.09 

H 

-1.X10' 10 

M 

1.X10' 3 

D 

-0.07 

I 

1.X10’ 6 

N 

2.X10’ 3 

E 

-0.05 

J 

5.X10' 5 




(b) Pressure Contour 


Label P 

Label P 

Label P 

A -900. 
B -650. 
C '400 . 

D -200. 

E 0. 

F 400. 

G 1000. 

H 3000. 


* i/>: Stream Function 
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The penalty method with the pressure interpolation polynomials of the form 
(l,x,y) yielded slightly distorted pressure contour lines near the top region of the 
cavity [Fig. 2(b)]. It was found that the distorted pressure contour lines had been 
caused by the ill-conditioned pressure matrix (M ) and the coefficients of the pres- 

sr 

sure interpolation polynomials. The entries of the pressure matrix were different by 
several orders of magnitude and the coefficients of the pressure interpolation poly- 
nomials were approximately 10 orders of magnitude different for the high aspect ratio 
fine grids located along the boundary of the cavity. However, these distorted pres- 
sure contour lines may disappear if a different pressure averaging technique and/or 
a different plotting package are used. 

The error norm (e.) versus number of iterations for each flow variable is shown 

in Figure 3. It can be seen in Figure 3 that the velocity-pressure integrated method 
and the penalty method with the pressure interpolation polynomials given as equations 
(15) and (17) yielded almost identical convergence rate for both velocity and pressure. 
But the penalty method with the pressure interpolation polynomials of the form (1 ,x,y) 
exhibited significantly degenerated convergence rate for pressure. It was found that 
the ill-conditioned pressure matrix (M^) and the computer round-off error were 

responsible for the degenerated convergence rate for pressure. For this case, the 
pressure does not seem to converge at all as the number of iterations were increased 
behond approximately 10 (Fig. 3). 



: Velocity- Pressure Integrated Method, 

: Penalty method with eq. (15), 

: Penalty Method with eq. (16). 

: Penalty method with eq. (17). 


Figure 3. Error norm versus number of iterations for cavity flow. 
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All the methods yielded almost identical velocity vectors, which can be seen 
from the stream line contours shown in Figures 1 and , 2. The horizontal velocity 
profiles at x = 0.5, obtained by using the velocity- pressure integrated method, com- 
pared favorably with those of Schreiber and Keller [18] and Ghia et al. [19] (see 
Reference 2). To solve the same cavity flow for Reynolds number of 10,000, 180 x 
180 grid points have been used in Schreiber and Keller [18]; and 129 x 129 grid 
points, in Ghia et al. [19]. 


3.2 Backward-Facing Step Flow 

A laminar backward- facing step flow is considered below. The experimental 
data can be found in Armaly et al. [22]. In the following discussion, the Reynolds 
number is based on the hydraulic diameter (D = 0.0104 m) and the bulk velocity (V = 
0.6667 m/sec) at the inlet. The experimental data showed that there exists only one 
recirculation zone at the down- stream region of the backward- facing step for the 
Reynolds number less than approximately 450. As the Reynolds number was increased 
beyond approximately 450, a second pressure driven recirculation zone appeared at 
the top wall of the channel. 

The Reynolds numbers considered in the present study were 410, 420, 430, 

440, and 500. These various Reynolds numbers have been obtained by varying the 
fluid viscosity. The velocity profile of a fully developed channel flow has been 
applied at the inlet boundary and the vanishing normal stress boundary condition has 
been prescribed at the exit boundary. The trivial solution (u = v = p = 0) has been 
used as an initial guess for all the Reynolds number cases considered. The quality 
of the solutions for all the methods remained unchanged after 50 iterations. A com- 
plete set of computational results obtained by using the velocity-pressure integrated 
method for the Reynolds number of 100 through 900 can be found in Reference 2. 

The finite difference computations of the same backward- facing step flow can be 
found in Armaly et al. [22] and Kim and Moin [23], among many others. 

All the methods considered herein yielded almost identical velocity vectors for 
Reynolds numbers of 430 and 440, respectively, and predicted the existence of the 
pressure driven recirculation zone at the top wall for Re >. 440. The streamline 
contours for Reynolds numbers of 430 and 440 are shown in Figure 4. The stream 
function (ip) has been normalized using a relationship given as [24]: 


f = <p/(U 


max 


h) 


( 22 ) 


where Y is the normalized stream function, U is the maximum velocity at the inlet, 

m ax 

and h is the step height. The streamline contour labels are given in Table 2. 

The streamline and pressure contours for Re = 500 obtained by using the 
penalty method with the pressure interpolation polynomials given in equations (15) 
and (16) are shown in Figures 5 and 6, respectively. The pressure (p) has been 
normalized using a relationship given as [ 24] : 

P = P /(pU max 2/2) (23) 
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(b) x/h = 38 

Figure 4. Streamline contours for backward- facing step flow, 

(a) Re = 430, and (b) Re = 440. 

TABLE 2. CONTOUR LABEL FOR BACKWARD-FACING STEP FLOW 


(a) Streamline Contour 


Label 


Label 


Label 


A 

-0.0408 

E 

0 . 0000 

I 

0.4082 

B 

-0.0305 

F 

0.0204 

J 

0.6122 

C 

-0.0102 

G 

0.1020 

K 

0.707483 

D 

-0.0020 

H 

0.2041 




(b) Pressure Contour 


Label 


-0.0027 

0.0017 

0.0111 

0.0314 

0.0755 



0.1333 

0.1556 

0.1778 

0.1955 

0.2011 


Label 


0 . 2044 
0.2058 
0.2069 











where P is the normalized pressure. The pressure contour labels are given in Table 
2. The reference pressure of p = 0 has been assigned at the concave corner of the 
step. The streamline contours obtained by using the velocity- pressure integrated 
method and the penalty method with the pressure interpolation polynomials given as 
equation (17) were identical to the streamline contour shown in Figure 5. The penalty 
method with the pressure interpolation polynomials of the form (l,x,y) yielded a 
severely distorted pressure contour (Fig. 6). 



CDE FGHIJK KJ I H G F 



Figure 5. Backward -facing step flow, penalty method with equation (15), 

(a) streamline, and (b) pressure. 
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Figure 6. Backward- facing step flow, penalty method with equation (16), 

(a) streamline, and (b) pressure. 
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The error norm versus number of iterations for each flow variable is shown in 
Figure 7. Both the velocity-pressure integrated method and the penalty method with 
the pressure interpolation polynomials given as equations (15) and (17) yielded rapidly 
convergent solutions . The pressure interpolation polynomials of the form (l,x,y) 
exhibited degenerated convergence rate for pressure, due to the same reasons as 
have been discussed in the previous cavity flow case. 



Figure 7. Error norm versus number of iterations for backward- facing 
step flow, notations are the same as in Figure 3. 



The present computational results compared favorably with the experimental 
data [22] as well as the finite difference computational results of Kim and Moin [ 23] ^ 
see Reference 2. In Kim and Moin [23], the exit boundary has been located at 30 
step heights downstream of the expansion corner and 101 x 101 grid points have been 
used. 


3.3 Flow Through a Nest of Cylinders 

Flows through a nest of cylinders can be found in a number of engineering 
applications such as the Space Shuttle Main Engine - Main Injector Assembly (SSME- 
4 MIA) and the heat exchangers in nuclear reactors (see Reference 25 for more details). 

However, these flows began to be solved numerically only very recently. These are 
a finite element computation of a two-dimensional laminar flow through a nest of 
■* cylinders [26] and a body-fitted grid finite difference computation of a three- 

dimensional laminar flows through a nest of cylinders [27]. Neither experimental 
data nor detailed computational results are available for these flows as yet. 

A laminar flow through a nest of cylinders at a Reynolds number of 40 is con- 
sidered below (Fig. 8). The Reynolds number is defined as Re = pUD/y, where 
U = 1 is the free stream velocity and D = 1 is the diameter of a cylinder. The inlet 
boundary has been located at three diameters upstream of the forward stagnation 
point of the first column of cylinders; and the exit boundary at 41 diameters down- 
stream of the inlet boundary. A uniform velocity profile has been used as the inlet 
boundary condition. The vanishing normal stress boundary condition has been pre- 
scribed at the exit boundary; and the symmetric boundary condition, at the top and 
the bottom of the computational domain. The computational domain has been discre- 
tized by 1024 quadratic elements with 4369 nodes. The finite element mesh in the 
vicinity of the nest of cylinders is shown in Figure 8. The trivial solution (u = v = 
p = 0) has been used as an initial guess. 



Figure 8. Flow through a nest of cylinders, grid in the 
vicinity of the nest of cylinders. 

The streamline and pressure contours obtained by using the velocity- pressure 
integrated method and the penalty method with the pressure interpolation polynomials 
given as equations (15) and (16) are shown in Figures 9, 10, and 11, respectively. 
The pressure has been normalized using a relationship given as P = p/(pu2/2), where 
U = 1 is the reference velocity at the inlet boundary. An arbitrary reference pres- 
sure (p = 0.0) has been assigned at the forward stagnation point of the first column 
of cylinders. The streamline contour label is given in Table 3. In Figures 9 through 
11, the minimum and maximum normalized pressures (P) are -20.0 and 0.0, respec- 
tively; and the incremental normalized pressure (AP) between the contour lines is 
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(a) 



(b) 


Figure 9. Flow through a nest of cylinders, velocity -pressure integrated 
method, (a) streamline, and (b) pressure. 





(b) 

Figure 10. Flow through a nest of cylinders, penalty method with 
equation (15), (a) streamline , and (b) pressure. 
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(a) 



(b) 


Figure 11. Flow through a nest of cylinders, penalty method with 
equation (16), (a) streamline, and (b) pressure. 

TABLE 3. STREAMLINE CONTOUR LABEL FOR FLOW THROUGH 

A NEST OF CYLINDERS 



equal to 1.0. The streamline and pressure contours obtained by using the penalty 
method with the pressure interpolation polynomials given as equation (17) was identi- 
cal to those shown in Figures 9 and 10. The penalty method with the pressure inter- 
polation polynomials given as equation (16) yielded severely distorted pressure contour 
lines for the same reasons as have been listed previously (Fig. 11). 

The error norm versus number of iterations for each flow variable is shown in 
Figure 12. The practically convergent solutions have been obtained after approxi- 
mately 15 iterations for all the cases. The velocity-pressure integrated method 
yielded uniformly convergent solution as before. The penalty method yielded rapidly 
convergent solutions as the velocity -pres sure integrated method at earlier iterations. 
For the arbitrary distorted quadrilateral elements with high aspect ratio, the adverse 
effect of the ill-conditioned pressure matrix and the computer round-off error became 
so severe that only the velocity-pressure integrated method yielded uniformly con- 
vergent pressure as the number of iterations was increased. 
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NO. OF ITERATIONS 

Figure 12. Error norm versus number of iterations for flow through 
a nest of cylinders, notations are the same as in Figure 3. 
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3.4 Channel Flow with an Internal Blockage 

There exists a controversy over the use of upwinding techniques for convec- 
tion dominated flows. Use of the upwinding techniques has been partly based on an 
argument that the discrete finite element system of equations become ill-conditioned 
as the grid Reynolds number is increased ; and it has been claimed that use of such 
upwinding techniques is the best approach to suppress the numerical wiggles when 
coarse grid is used [28]. The opponents claimed that the numerical diffusion intro- 
duced by use of such upwinding techniques may obscure the physical diffusion 
process and the computational results may not be accurate [ 29] . 

A channel flow with an internal blockage is considered below to further investi- 
gate the cause of numerical wiggles. The finite element mesh for the full computa- 
tional domain is shown in Figure 13. The velocity profile of a fully developed channel 
flow has been used as the inlet boundary condition; and the vanishing normal stress 
has been prescribed at the exit boundary. The trivial solution (u = v = p = 0) has 
been used as an initial guess. The streamline and pressure contours obtained by 
using the velocity- pressure integrated method and the penalty method with the pres- 
sure interpolation polynomials given as equations (15) and (16) are shown in Figures 
14, 15, and 16, respectively. The convergence history for each of the flow variables 
is given in Figure 17. Again, only the velocity- pressure integrated method yielded 
uniformly convergent solution as the number of iterations was increased. 



Figure 13. Computational domain and finite element grid for 
channel flow with an internal blockage . 




(b) 

Figure 14. Channel flow with an internal blockage, velocity-pressure 
integrated method, (a) streamline, and (b) pressure. 
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(b) 

Figure 15. Channel flow with an internal blockage, penalty method with 
equation (15), (a) streamline, and (b) pressure. 




Figure 16. Channel flow with an internal blockage, penalty method with 
equation (16), (a) streamline , and (b) pressure. 







It can be seen in Figures 14 through 16 that there exists a steep pressure 
gradient in the forward corner region of the blockage. The same flow has been 
solved using four different, locally refined, coarse grids (Fig. 18). It can be seen 
in Figure 19(b) that the grid refinement in the local high Reynolds number region 
was not helpful to suppress the numerical wiggles. On the other hand, any grid 
refinement in the steep pressure gradient region suppressed the numerical wiggles 
significantly [Fig. 19(c)-(d)] . For this flow case, it can be concluded that the 
numerical wiggles have been caused by the coarse grids [Fig. 18(a)-(b), which could 
not resolve the steep pressure gradient in the forward corner region of the blockage. 
Note that the local Reynolds number in the high pressure gradient region is suffi- 
ciently small compared with that of the upstream region. Use of an upwinding tech- 
nique has been partly justified based on the assumption that the high grid Reynolds 
number is responsible for numerical wiggles. However, these computational results 
suggest that the high grid Reynolds number was less responsible for the numerical 
wiggles than the steep gradient of a flow variable, which turned out to be the pres- 
sure for this flow case, in the forward corner region of the blockage. 



(a) 



(b) 



(c) 
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Figure 18. Local grid refinement for channel flow 
with an internal blockge. 
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IV. CONCLUSIONS AND DISCUSSION 


A comparative study of the velocity- pressure integrated finite element method 
and the consistent penalty finite element method has been presented. The penalty 
method with the pressure interpolation polynomials given as equations (15) and (17) 
yielded uniformly convergent solutions. The convergence rate was equal to that of 
the velocity-pressure integrated method. The penalty method with the pressure 
interpolation polynomials of the form (l,x,y) exhibited slightly degenerated conver- 
gence rate. 

It was found that all of the methods yielded almost identical computational 
results for the velocity. However, the pressure interpolation polynomials of the form 
(l,x,y) yielded severely distorted pressure contours for the example flow cases. 

The distorted pressure contours had been caused by the ill-conditioned matrix of 
the discrete penalized conservation of mass equation. The penalty methods required 
slightly smaller computational time than the velocity-pressure integrated method. 
: HoweverT-t : he-dif-ferenee-was-Thsighi : ficahf.-^T-h'e-velocity^pressure-integrated-;method 
would be preferable over the penalty methods, for its uniform convergence behavior 
for pressure. 

It has been shown that any of the finite element methods considered in this 
report could capture the subtle pressure driven recirculation zones for high Reynolds 
number flows. The computational results compared favorably with experimental data 
and/or fine grid finite difference computational results. 

For the example problems considered herein, a relatively small number of grid 
points, compared with the fine grid finite difference computations of the same example 
flows, were required to resolve the details of the flow field. It was found that no 
upwinding technique was necessary to obtain computational results which were free 
of numerical wiggles for high Reynolds number flows. 
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APPENDIX I 

FINITE ELEMENT COMPUTER PROGRAM (NSFLOW/P) FOR 
INCOMPRESSIBLE, LAMINAR FLOWS 


PRECEDING PAGE BLANK NOT FILMED 
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o o o o noo non 


********************************* TOP OF DATA ********************************* 
C 

C********l*********2*********3*********4*********5*********6*** 

C PROGRAM MAIN 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

CKARACTER*4 TITLE, IWORD 
DIMENSION TITLE(15) , IWORD(IO) 

DATA IWORD / 'INIT', 'PREP', '****', 'PROC', 'CONT', 

‘ **** ' f ' **** ' f ' **** ' ) ' **** ' t < end ' / 

DATA MAXNOD , MAXELM , MAXDOF , MXFRON /4227, 1027, 11527, 167/ 

C 

101 CONTINUE 

READ (5, 501) TITLE 
WRITE(6 , 601) TITLE 
501 FORMAT (20A4) 

601 FORMAT(/2X, 20A4) 

C 

102 CONTINUE 

READ (5, 501) TITLE 
WRITE (6, 601) TITLE 
DO 103 K-1,10 

IF(TITLE(2) . EQ. IWORD (K) ) GO TO 105 

103 CONTINUE 

104 WRITE(6 , 602) 

602 FORMAT (2X, 'TERMINATED IN SUB-PREP FOR INPUT DATA ERROR') 

STOP 

C 

105 CONTINUE 

GO TO (1,2, 3, 4,5, 6,7,8,9,10), K 

INITIALIZE DIMENSIONED VARIABLES 

1 CONTINUE 

CALL INITAL(MAXNOD , MAXELM , MAXDOF , MXFRON) 

GO TO 102 

PREPARE INPUT DATA 

2 CONTINUE 

CALL PREP (MAXNOD, MAXELM, MAXDOF, MXFRON) 

GO TO 102 

3 CONTINUE 
GO TO 104 

UNSTEADY FLOW SOLVER 

4 CONTINUE 

CALL PROCES (MAXNOD, MAXELM, MAXDOF, MXFRON) 

GO TO 101 

5 CONTINUE 
GO TO 101 
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6 CONTINUE 

7 CONTINUE 

8 CONTINUE 

9 CONTINUE 
10 CONTINUE 

STOP 

END 

C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE INITAL(MAXNOD , MAXELM , MAXDOF , MXFRON ) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CPRS/ PELEM(4, 1027) , PBCDAT, IPNOD(2) , IPDOF 
COMMON /CGRID/ X(4227 , 3) ,NODES(27 , 1027) 

COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 

C 

DO 10 KDIM=1 , 3 
DO 10 KNODE=l , MAXNOD 
X(KNODE.KDIM) =0. 

10 CONTINUE 
C 

DO 30 KELEM-1, MAXELM 
DO 30 KPE=1 , 27 
NODES (KPE , KELEM) =0 
30 CONTINUE 
C 

DO 50 KPROB=l ,10 
DO 50 KNODE=l , MAXNOD 
I BCA( KNODE, KPROB) = 0 
ADBC ( KNODE , KPROB ) = 0. 

A ( KNODE , KPROB ) = 0. 

50 CONTINUE 
C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6 
BLOCK DATA BLKDAT 
C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE , NELEM , NPE , NPRE , NDIM , NEDOF , IFLOW , 

IAXSY , IELF 

COMMON /CGAUL/ CLXKS (4,4), CLW (4,4), NGAUS 
COMMON /CGAUS/ EXKS ( 3 , 64) , EW( 64) , MGAUS 
COMMON /CINDX/ INDXF(27 , 3 , 15) , INDXP(27 , 15) 

COMMON /CITER/ CNVCF(IO) , ERROF(IO) , RELAX( 10) , ITERE , MAX1T 
COMMON /CLSCF/ XKSNOD(27 , 3 , 11) ,TM(4 , 4) 

COMMON /CMATE/ BFX(3) ,DENSY,VISCY, PECLET 
COMMON /CPNLT/ CXKPN(4 ,3,15), CWPN(4 , 15) , PTNUM 
COMMON /CPROB/ IA(10),IPL0T 
C 

DATA CLXKS(l.l) /0./, 

CLW (1,1) /2 . /> 

(CLXKS (K, 2) ,K=1 , 2) /-0 . 5773502692 , 0.5773502692/, 
(CLW(K, 2) ,K=1 , 2) / 1., 1./, 

(CLXKS (K , 3 ) , K=1 , 3) /-0 . 7745966692 , 0., 0.7745966692/, 



o o 


C 


C 


(CLW(K, 3) ,K=l,3)/0. 5555555556, 0.8888888889, O'. 5555555556/, 
(CLXKS(K,4) ,K=1,4) /-0 . 8611363116 , -0.3399810436, 

0.3399810436, 0.8611363116/, 

(CLW(K,4) ,K=1 ,4) / 0.3478548451, 0.6521451549, 

0.6521451549, 0.3478548451/ 

DATA (CXKPN(K, 1,4) ,K=1 , 3) /0., -0.70710678, 0.70710678/, 

(CXKPN(K, 2 ,4) , K=1 , 3 ) /O . 81649658 , -0 .40824829 , -0 . 40824829/ , 
(CWPN(K,4) ,K=1, 3) /l. 33333333, 1.33333333, 1.33333333/ 

DATA (XKSNOD(K, 1,6) ,K=1 , 4) / -1., 1., 1., -1./, 

(XKSN0D(K, 2,6), K=1 , 4) / -1., -1., 1., 1./, 

(XKSN0D(K, 1,8), K=1 ,9) / -1 . , 0 . , 1 . , 1 . , 1 . , 0 . , - 1 . , -1 . , 0 . / , 
(XKSNOD(K, 2 , 8) , K— 1,9) / -1 . , -1 . , -1 . , 0 . , 1 . , 1 . , 1 . , 0 . , 0 . / 


--- IFL0W=5 FOR 2-D CONSISTENT PENALTY METHOD 

DATA ( INDXF (KPE ,1,5) , KPE=1 , 9 ) / 1, 3, 5, 7, 9,11,13,15,17/, 

( INDXF (KPE ,2,5) , KPE=1 , 9 ) / 2, 4, 6, 8,10,12,14,16,18/ 

C --- IFL0W=6 FOR 2-D IMPROVED CONSISTENT PENALTY METHOD 

DATA ( INDXF (KPE ,1,6) , KPE=1 , 9 ) / 1, 3, 5, 7, 9,11,13,15,17/, 

( INDXF ( KPE , 2 , 6 ) , KPE=1 , 9 ) / 2, 4, 6, 8,10,12,14,16,18/ 

C 

DATA ((TM(I,J) ,J=1,4) ,1=1,4) / 0.25, 0.25, 0.25, 0.25, 

-0.433012701,-0.433012701, 0.433012701, 0.433012701, 
-0.433012701, 0.433012701,-0.433012701, 0.433012701, 

0.75, -0.75, -0.75, 0.75/ 

C 

DATA I FLOW , IAXSY , I PLOT /3*0/, 

NPRE , MGAUS , NGAUS , MAXIT /4*0/ 

DATA VISCY , DENSY , PECLET /3*0./ 

DATA XKAPA , EWALL , CMUF , TKRWAL , TKREXT , XMLTH /6*0./ 

DATA ( IA(K) , K=1 , 10) /10*0/ 

DATA ERROF/10*0 . / 

END 

C 

C********i*********2*********3*********4*********5*********6*** 
SUBROUTINE DATLIB 
C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE , NELEM , NPE , NPRE , NDIM , NEDOF , IFLOW , 

IAXSY, IELF 

COMMON /CGAUL/ CLXKS(4,4) ,CLW(4,4) , NGAUS 
COMMON /CGAUS/ EXKS (3 , 64) , EW(64) , MGAUS 
COMMON /CPNLT/ CXKPN(4 , 3,15) , CWPN(4 , 15) , PTNUM 
COMMON /CPROB/ IA(10),IPLOT 

DIMENSION LIBELF(15) ,LIBELH(6) .LIBNPE(ll) ,LBNPRE(15) 

C 

DATA (LIBELF(IFL) , IFL=1 , 15) /0, 0,0, 0,8, 8, 0,0, 0,0, 0,0, 0,0,0/, 
(LIBNPE(IEL) , IEL=1 , 11) /0, 0,0, 0,0, 4, 8, 9, 0,0, 0/, 
(LBNPRE(IFL) , IFL=1,15) /0, 0,0, 0,3, 3, 0,0, 0,0, 0,0, 0,0,0/ 

C 

IELF=LIBELF( IFLOW) 

NPRE=LBNPRE (IFLOW) 

NPE. =LIBNPE(IELF) 

C 
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LGAUS=0 

MGAUS=NGAUS**NDIM 
GO TO (10,20,30) , NDIM 
10 WRITE(6 , 601) NDIM 
STOP 

601 FORMAT (2X, 'TERMINATED AT SUB-DATLIB NDIM=',I2) 

C 

20 CONTINUE 

DO 2 1=1 , NGAUS 
DO 2 J=l, NGAUS 
LGAUS=LGAUS+1 

EXKS ( 1 , LGAUS )=CLXKS ( I , NGAUS ) 

EXKS ( 2 , LGAUS ) =CLXKS ( J , NGAUS ) 

EW(LGAUS) =CLW(I , NGAUS )*CLW(J , NGAUS) 

2 CONTINUE 
GO TO 100 

C 

30 CONTINUE 

DO 3 1=1, NGAUS 
DO 3 J=l, NGAUS 
DO 3 K=l, NGAUS 
LGAUS=LGAUS+1 

EXKS ( 1 , LGAUS ) =CLXKS ( I , NGAUS ) 

EXKS ( 2 , LGAUS ) =CLXKS ( J , NGAUS ) 

EXKS ( 3 , LGAUS ) =CLXKS (K , NGAUS ) 

EW( LGAUS ) =CLW( I , NGAUS ) *CLW ( J , NGAUS )*CLW(K, NGAUS ) 

3 CONTINUE 
C 

100 CONTINUE 
C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE PREP (MAXNOD , MAXELM , MAXDOF , MXFRON) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

CHARACTER*4 TITLE , ICNTRL, IWORD 

COMMON /CDESC/ NNODE , NELEM , NPE , NPRE , NDIM , NEDOF , I FLOW , 
IAXSY, IELF 

COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , 1BCA(4227 , 10) 
COMMON /CFRON/ MFRONF 

COMMON /CGAUL/ CLXKS (4 , 4) , CLW(4 , 4) , NGAUS 
COMMON /CGAUS/ EXKS (3 , 64) , EW( 64) ,MGAUS 
COMMON /CGRID/ X(4227 ,3) ,NODES(27, 1027) 

COMMON /CINDX/ INDXF(27 , 3 , 15) , INDXP(27 , 15) 

COMMON /CITER/ CNVCF(IO) ,ERROF(10) .RELAX (10) ,ITERE,MAXIT 
COMMON /CMATE/ BFX(3) .DENSY.VISCY, PECLET 
COMMON /CPNLT/ CXKPN(4 , 3,15), CWPN(4, 15) , FTNUM 
COMMON /CPROB/ IA(10),IPL0T 

COMMON /CPRS/ PELEM(4 , 1027) , PBCDAT, IPNOD(2) , IPDOF 
COMMON /CWIND/ WHI (27) , DWHX(27 , 3) 

DIMENSION TITLE(15),IWORD(30) 

DATA IWORD/ 'DESC' , ' CNTL' , 'ELEM', 'NODE', 'MATE', 
'****' f '****' t 'ITER' , '■****' t '***•*' j 



- 

' IA01 ' , 

' IA02 ' , 

' IA03 ' , 

' IA04 ' , 

» **** < 

» 

- 

' **** ' 

1 

' **** ' 

f 

' **** ' t 

' **** ' 

y 

'****» 

* 

- 

'****' t 

'****' 

1 - 

'****< t 

• **** • 

f 

' **** ' ( 

- 

' **** ' 

> 

' INCL' , 

• **** ' 

» 

1 **** • 

y 

'END '/ 


101 CONTINUE 

READ(5 , 501) ICNTRL, TITLE 
WRITE (6, 601) ICNTRL, TITLE 
501 FORMAT (20A4) 

601 FORMAT (/2X,20A4) 

C 

DO 102 K=1 ,30 

I F (ICNTRL . EQ . IWORD ( K ) ) GO TO 105 

102 CONTINUE 

103 WRITE (6 , 602) 

WRITE(6 , 601) ICNTRL, TITLE 

602 FORMAT (2X, 'TERMINATED IN SUB-PREP FOR INPUT DATA ERROR') 
STOP 


105 CONTINUE 

GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, 
21,22,23,24,25,26,27,28,29,30), K 


C 


1 CONTINUE 

READ (5, 501) TITLE 
WRITE(6 , 603) TITLE 

READ. (5',*) I FLOW , NDIM , NGAUS , MFRONF 

WRITE(6 , 605) I FLOW , ND IM , NGAUS , MFRONF 
IF (MFRONF. GT.MXFRON) GO TO 103 
603 FORMAT (2X.20A4) 

605 FORMAT (4X, 'IFLOW=' ,12, 2X, 'NDIM =' , 12 , 

2X, 'MFRONF=' ,14) 


2X, 'NGAUS=' ,12, 


CALL DATLIB 

WRITE (6, 606) I ELF , NPE , NPRE , MGAUS 
606 FORMAT(4X, 'IELF=' ,12, 2X, 'NPE-' ,12 , 2X , ' NPRE= ',12, 

2X, 'MGAUS=' ,12) 

GO TO 101 

2 CONTINUE 

READ(5 , 501) TITLE 
WRITE(6, 603) TITLE 
READ (5,*) NNODE , NELEM , IAXSY , I PLOT 
WRITE (6, 6 10) NNODE, NELEM, IAXSY, IPLOT 
I F ( NNODE . GT . MAXNOD . OR . NELEM . GT . MAXELM ) THEN 
WRITE(6 , 611) NNODE, NELEM, MAXNOD, MAXELM 
STOP 
END IF 

610 FORMAT(2X, 'NNODE=' ,15, 2X, 'NELEM=' , 15 , 2X, ' IAXSY=' , 12, 

2X, 'IPLOT=' ,12) 

611 FORMAT (2X, 'TERMINATED IN SUB-PREP FOR NNODE=',I6, 

2X, 'NELEM=' ,16, 2X , ' MAXNOD= ' ,16, 2X, 'MAXELM=' , 16) 

GO TO 101 
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3 CONTINUE 

CALL RELEM( NODES , NELEM , NPE , MAXELM) 

GO TO 101 
C 

4 CONTINUE 

CALL RNODE(X,NNODE,NPE, IELF.NDIM, MAXNOD) 

GO TO 101 
C 

5 CONTINUE 

READ (5 , 501) TITLE 
WRITE(6 , 603) TITLE 

READ (5,*) VISCY, DENSY, (BFX(K) ,K=1,NDIM) 

WRITE(6 , 613) VISCY, DENSY, (BFX(K) ,K=1,3) 

613 FORMAT ( 2X, 'VISCY=' ,E12. 4, 2X, ' DENSY= ' , E12 .4 , 

/2X, 'BFX=' .3E12.4) 

GO TO 101 
C 

6 CONTINUE 
GO TO 101 

C 

7 CONTINUE 
GO TO 101 

C 

8 CONTINUE 

READ (5, 501) TITLE 
WRITE(6 , 603) TITLE 

READ(5 ,*) MAXIT, (RELAX(K) , K=1 , 10) , (CNVCF(K) , K-l , 10) 
WRITE(6 , 626) MAXIT, (RELAX(K) , K-l, 10) , (CNVCF(K) , K-l, 10) 
626 FORMAT ( 2X, 'MAXIT-' ,15, 

/2X, 'RELAX=' .5E12.4, /8X.5E12.4, 

/2X , ' CNVCF= ' , 5E12 . 4 , /8X.5E12.4) 

GO TO 101 
C 

9 CONTINUE 
GO TO 101 

C 

10 CONTINUE 
GO TO 101 

C 

11 CONTINUE 

CALL RINIT(A(1 , 1) .NNODE , MAXNOD ) 

CALL RBC1(ADBC(1 , 1) ,IBCA(1,1) , MAXNOD , NNODE) 

GO TO 101 
C 

12 CONTINUE 

CALL RINIT(A(1, 2) , NNODE, MAXNOD) 

CALL RBC1 (ADBC (1,2), IBCA(1,2) .MAXNOD , NNODE) 

GO TO 101 
C 

13 CONTINUE 

CALL RINIT(A(1, 3) , NNODE, MAXNOD) 

CALL RBC1( ADBC (1,3) ,IBCA(1,3) .MAXNOD, NNODE) 

GO TO 101 
C 
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C 1 2 3 4 --5 6- 

14 CONTINUE 

READ (5 ,*) PBCDAT, (IPNOD(K) ,K=1 , 2) 

WRITE(6 , 635) PBCDAT, (IPNOD(K) ,K=1,2) 

635 FORMAT (2X, 'PRESSURE B.C. DATA PBCDAT=' ,E12 .4, 

2X, ' IPNOD(l-2)=' ,217) 

GO TO 101 

C 1 2 3 4 5 6 

C 

15 CONTINUE 
GO TO 101 

C 

C- -1 2 3 --4 5 --6 

16 CONTINUE 
GO TO 101 

C 

C ---1 2 3 7 5 6 

17 CONTINUE 
GO TO 101 

C 

18 CONTINUE 
GO TO 101 

C 

19 CONTINUE 
GO TO 101 

C 

20 CONTINUE 
GO TO 101 

C 

21 CONTINUE 
GO TO 101 

C 

22 CONTINUE 
GO TO 101 

C 

C TURBULENCE DATA 
C 

23 CONTINUE 
GO TO 101 

C 

24 CONTINUE 
GO TO 101 

C 

25 CONTINUE 
GO TO 101 

C 

26 CONTINUE 
GO TO 101 

C 

C INCLUDE RE-START DATA 

C 

27 CONTINUE 

CALL FEMDAT ( A , ADBC , X , PBCDAT , NODES , IBCA , IPNOD , NPE , NNODE , 
NELEM , MAXNOD , MAXELM) 
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GO TO 101 


1 

28 CONTINUE 
GO TO 103 

i 

--2 

_ _ 2 

3- 

3 

4 

3 

6 

^ 

29 CONTINUE 
RETURN 

--1 

-.2 

3 


5 

6--- 


30 CONTINUE 
C 

PTNUM= (VISCY/DENSY) *1 . E+10 
C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE RNODS (X , NNODE , NPE , I ELF , NDIM , MAXNOD) 

C-X- IMPLICIT REAL *8 (A-H.O-Z) 

CHARACTERS TITLE 

DIMENSION X (MAXNOD , 3) , DELX(197 , 3) ,XNOD(27,3) ,NKS(3) ,CXKS(3) , 
PHI (27) , DPHI( 27 , 3) , WHI(27) , DWHI (27 , 3) .TITLE (15) 

C 

READ (5, 501) TITLE 
WRITE(6 , 601) TITLE 
501 FORMAT (20A4) 

601 FORMAT (2X.15A4) 

C 

READ (5,*) NBLOC 
WRITE(6 , 605) NBLOC 

605 FORMAT(2X, 'SUB-RNODE NBLOC=',I2) 

C 

DO 7000 IBLOC-1, NBLOC 
READ(5 , 501) TITLE 
WRITE(6 , 601) TITLE 
READ (5 , *) METHOD 
WRITE (6, 606) METHOD 

.606 FORMAT (4X, 'GRID GENERATION METHOD=',I3) 

GO TO (1000,2000) METHOD 
C 

1000 CONTINUE 

READ (5 , 501) TITLE 
WRITE(6 , 601) TITLE 
READ ( 5 , * ) NODG1, INCRX, INCRY, INCRZ 
WRITE (6, 640) NODG1 , INCRX, INCRY, INCRZ 
640 F0RMAT(4X, 15 , 2X,I3, 2X,I3, 2X,I3) 

READ(5 , 501) TITLE 
WRITE(6 , 601) TITLE 
DO 10 KDIM=1,3 

READ ( 5 , * ) NDAT , (DELX( IKE , KDIM) , IKE=1 ,NDAT) 

WRITE(6 , 642) NDAT, (DELX( IKE, KDIM) , IKE=1 ,NDAT) 

NKS (KDIM) =NDAT 
IF (NDAT . GT . 197) THEN 
WRITE(6 , 645) 

STOP 
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END IF 

10 CONTINUE 

642 FORMAT (2X„'NDAT=' ,15, 20(/4X, 5F10 . 7) ) 

645 FORMAT (2X, 'INPUT DATA ERROR FOR NDAT IN SUB-RNODE') 
LLINE=NKS ( 1 ) 

MLINE=NKS ( 2 ) 

NLINE=NKS(.3) 

DO 15 KLINE=1 , NLINE 
DO 15 JLINE=1 .MLINE 
DO 15 ILINE=1 , LLINE 

KNODE=NODGl+ (ILINE - 1 ) *INCRX+ (JLINE - 1 ) *INCRY+ ( KLINE - 1 ) *INCRZ 
X(KNODE, 1)=DELX (ILINE, 1) 

X(KNODE, 2 )=DELX( JLINE , 2 ) 

X(KNODE , 3 ) =DELX ( KLINE , 3 ) 

15 CONTINUE 
GO TO- 7000 

2000 CONTINUE 

READ (5, 501) TITLE 
WRITE (6, 601) TITLE 

READ (5,*) ( (XNOD (KPE , KDIM) ,KDIM=1 , NDIM) ,KPE=1,NPE) 

DO! 22 KPE=1 ,NPE 

WRITE(6 , 610) KPE, (XNOD (KPE, KDIM) , KDIM=1 ,NDIM) 

22 CONTINUE 

610 FORMAT ( 4X , ' KPE= ' ,12, 2X.3E12.4) 

READ(5 , 501) TITLE 
WRITE (6 ,601).' TITLE 
READ (5,*) NODG1 , INCRX , INCRY , INCRZ 
WRITE ( 6,612 ) NODG1 , INCRX .INCRY , INCRZ 
612 FORMAT (4X, 15 , 2X.I3, 2X,I3, 2X,I3) 

READ(5 ,501) TITLE 
WRITE(6 ,601) TITLE 
DO 25 KDIM=1 , 3 

READ (.5,*) NDAT, ( DELX( IKE, KDIM) , IKE=1 , NDAT) 

WRITE (6, 6 14) NDAT, (DELX( IKE, KDIM) , IKE-1, NDAT) 

NKS (KDIM) =NDAT 
IF (NDAT . GT . 197) THEN 
WRITE (6, 6 20) 

STOP 

ENDIF 

25 CONTINUE 

614 FORMAT ( 2X , ' NDAT= ',15, 10(/4X, 10F5 . 2) ) 

620 FORMAT (2X, 'INPUT DATA ERROR FOR NDAT IN. SUB-RNODE ' ) 
LLINE=NKS (1) 

MLINE=NKS(2) 

NLINE=NKS (3) 

DO 45 KLINE=1, NLINE 
DO 45 JLINE=1 , MLINE 
DO 45 ILINE=1, LLINE 
CXKS ( 1 ) =DELX ( ILINE , 1 ) 

CXKS (2)=DELX(JLINE , 2) 

CXKS (3 )=DELX (KLINE, 3) 
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GOTO (31,31,31,31,31, 31,31,38,31,31, 31), I ELF 

31 CONTINUE 

WRITE(6 ,630) I ELF 

630 FORMAT ( 2X, ' SUB- RNODE IELF=',I2) 

STOP 

C 

38 CALL SHAP23(PHI ,DPHI , CXKS ,NPE) 

C 

42 CONTINUE 

KNODE=NODGl+ ( ILINE-1 ) *INCRX+ ( JLINE - 1 ) *INCRY+(KLINE - 1 ) *INCRZ 

DO 44 KDIM=1 , NDIM 

X(KNODE,KDIM)=0. 

DO 44 KPE=1 ,NPE 

X(KNODE,KDIM)=X(KNODE,KDIM)+XNOD(KPE,KDIM)*PHI(KPE) 

44 CONTINUE 

45 CONTINUE 


7000 CONTINUE 
RETURN 
END 
C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE RELEM (NODES , NELEM , NPE , MAXELM) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

CHARACTERS TITLE 

DIMENSION NODES (2 7, MAXELM) , NEL( 3 ) , INCREL( 3 ) , INCNOD (27,3) , 
TITLE (15) 

C 

DO 1 KELEM=1 , NELEM 
DO 1 KPE==1 ,NPE 
NODES (KPE,KELEM)=0 
1 CONTINUE 
C 

READ (5, 501) TITLE 
WRITE(6 , 601) TITLE 
501 FORMAT ( 2 0A4) 

601 FORMAT ( 2X,1 5 A4) 

READ (5,*) NBLOC 
WRITE (6, 610) NBLOC 
610 FORMAT (4X, 'NBLOC=' ,16) 

C 

DO 100 IBL0C=1, NBLOC 
READ (5 , 501) TITLE 
WRITE(6 , 601) TITLE 

READ(5 ,*) IEL1 , (NODES (I PE , I ELI) , IPE=1 ,NPE) 

WRITE(6 , 620) IEL1 , (NODES(IPE , IEL1) , IPE=1 ,NPE) 

620 F0RMAT(4X, ' IEL1=' , 16 , 2X, ' NODES (KPE , IEL1)=' , 816 , 

3(/15X, ' NODES (KPE, IEL1)=' ,816)) 

READ (5, 501) TITLE 
WRITE(6 , 601) TITLE 
DO 10 KDIM=1 , 3 

READ ( 5 , * ) NEL(KDIM) , INCREL(KDIM) , ( INCNOD (K.KDIM) ,K=1,NPE) 
10 CONTINUE 
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NELX=NEL(1) 

NELY=NEL(2): 

NELZ=NEL(3) 

C 

DO- 50 IELZ=li , NELZ- 
DO 50 IELY=1 ,NELY 
DO 50 IEEX=“I.,.NELX 

KELEM=IEL1+(IELX-1)*INCREL(,1) + (IELY- I)*INCREL(2) 

+(IELZ- 1)*INGREL( 3) 

DO 30 KPE=1 , NPE. 

NODES (KPE , KELEM) =NODES (KPE , I ELI ) + ( I ELX - 1 ) *INCNOD ( KPE , 1 ) 
+ ( I ELY- 1 ) *INCNOD (KPE ,2) + ( IELZ -1.) *INCNOD (KPE', .3 ), 

30- CONTINUE 
50 CONTINUE, 

C. 

loo continue; 
c 

RETURN, 

END- 


C***.*****l**K******2********.*3*tf*******4***-****tf*5*********6*** 
SUBROUTINE. EINIT (AINIT, NNODE , MAXNOD) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

DIMENS ION: AINIT (MAXNOD )> ,,TTTLE( 1 5 ) 

C' 

READ (5, 501) TITLE’ 

WRITE (6., 601). TITLE 
501. FORMAT (I5A4) 

601 FORMAT (/2X,15A4). 

READ (-5.*) NREC 
WRITE ( 6 „ 6 10 )) NREC: 

IE (NREC . LE‘..0 )> RETURN 
610 FORMAT ( 2X , 'SUB- RTNIT NREC=* ,15) 

C. 

DO 20 IREC-1 , ; NREC 

READ (5,*) N1 „N2 , INCNOD-, ADATA. 

WRITE(6 , 620) N1 ,.N2 .INCNOD .ADATA 
620 FORMAT’( 5X ,. ' Nl= ' ,16 , 5X, 'N2=",I6, 5X„'TNCNOD=' ,I6 „ 

5X:„'ADATA-',E12..4) 

DO' 10- N=N1,.,N2.,, INCNOD 
AINIT (N)=ADATA 
10 CONTINUE 
20 CONTINUE. 

C- 

RETURN, 

end; 


* ***** ** 1 'k-li ii'k * * ** *2 * * * * *3 *•** ****** 4 ******* -irtc 5 * ******** 6 * * * 

SUBROUTINE. RBC1 ( DBCH , LDBC', MAXNOD-, NNODE ) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

GHARACTER*4 TITLE’, 

DIMENSION DBCH (MAXNOD) ,. LDBC (MAXNOD) ,TITLE(15) 

C 

DO! 10 KNODE=l,. NNODE 
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LDBC ( KNODE ) =0 
DBCH(KNODE)=0. 

10 CONTINUE 
C 

DBC DATA 

READ(5 , 501) TITLE 
WRITE(6 , 601) TITLE 
READ (5,*) NREC 
WRITE (6, 602) NREC 
IF(NREC.EQ.O) RETURN 
C 

WRITE(6 , 603) 

DO 40 IREC=1 , NREC 
READ(5 ,*) N1 , N2 , INCR , DUM 
WRITE (6, 604) N1 , N2 , INCR , DUM 
DO 30 K=N1 , N2 , INCR 
LDBC(K)=1 
DBCH(K)=DUM 
30 CONTINUE 
40 CONTINUE 
C 

WRITE(6 , 607) 

DO 60 KNODE=l , NNODE 

IF (LDBC (KNODE) .NE.O) WRITE(6,605) KNODE, LDBC ( KNODE ) , 

DBCH(KNODE) 

60 CONTINUE 
C 

RETURN 

501 FORMAT (20A4) 

601 FORMAT (/2X.20A4) 

602 FORMAT (5X, 'NO. OF INPUT DATA RECORD FOR DBC, NREC=',I5) 

603 FORMAT ( 5X , ' Nl-NODE N2-NODE INCREMENT DBC - DATA ' ) 

604 FORMAT (5X, 15 , 5X, 15 , 5X, 15 , 5X, Ell .4) 

605 FORMAT ( 5X , ' NODE= ' , 15 , 5X, ' LDBC=' ,I3,5X, ' DATDBC= ' ,E11.4) 

607 FORMAT (/2X, 'LIST OF D.B.C. DATA FROM SUB-RBC1') 

END 

C 

C********i*********2*********3*********4*********5*********6*** 
SUBROUTINE FEMDAT ( A , ADBC , X , PBCDAT , NODES , IBCA , IPNOD , NPE , 

NNODE , NELEM , MAXNOD , MAXELM) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

CHARACTER*4 TITLE 

DIMENSION A(MAXN0D , 10) , ADBC (MAXNOD, 10) ,X(MAXN0D,3) , 

NODES (27, MAXELM) , IBCA(MAXNOD , 10) ,IPNOD(2) ,TITLE(15) 
C 

READ(4 , 501) TITLE 
DO 50 KNOD-1, NNODE 

READ(4,*) KNODE, (X(KNODE.KDUM) ,KDUM=1,3) 

50 CONTINUE 
501 FORMAT (15A4) 

C 

READ(4 , 501) TITLE 
DO 52 KEL=1, NELEM 



READ(4,*) KELEM, (NODES (KPE,KELEM) ,KPE=1,NPE) 

52 CONTINUE 
C 

READ(4 , 501) TITLE 
DO 54 KNOD=I , NNODE 

READ(4 , *) KNODE, (IBCA(KNODE.K) ,K=1,3) 

54 CONTINUE 

READ(4 , 501) TITLE 

READ(4 , *) PBCDAT, (IPNOD(K) ,K=1,2) 

C 

READ(4 , 501) TITLE 

DO 55 KNOD=l , NNODE 

READ(4,*) KNODE, (A(KNODE,K) ,K=1,4) 

55 CONTINUE 

C 

READ(4 , 501) TITLE 
DO 57 KNOD=l, NNODE 

READ(4, *) KNODE, (ADBC (KNODE ,K) ,K=1 , 3) 

57 CONTINUE 
C 

DO 60 KNODE=l , NNODE 
A(KNODE , 4)=0 . 

60 CONTINUE 
C 

RETURN 

END 

C 

C********i*********2*********3*********4*********5*********6 
SUBROUTINE ISOPEL( IFLOW , IELF, IAXSY , NPE , NPRE , NDIM) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CELEM/ EX(27 , 3 ) , EA(27 , 10) , ENUT(64) , ENUT1 (8 ) , 
NODEL(27) , IELEM 

COMMON /CWIND/ WHI ( 2 7 ) , DWHX ( 2 7 , 3 ) 

COMMON /ESHAP/ PHI (27 ) , DPHX(27 , 3) , PSI ( 27 ) , DXDK(3 , 3 ) , 

DKDX (3,3) , CXKS ( 3 ) , DETJB , RADUS , IGAUS , JGAUS , LGAUS 
DIMENSION DPHI (27 , 3) ,DPSI(27,3) ,XGS(3) 

C 

GO TO (1,1, 1,1,1, 1,1, 8, 1,1, 1), IELF 
1 CONTINUE 

WRITE(6 , 608) IFLOW, IELF 

608 FORMAT ( 2X, ' SUB- I SOPEL IFLOW=',I2, 2X, ' IELF=' , 12) 

STOP 

C 

8 CALL SHAP23 (PHI, DPHI, CXKS, NPE) 

C 

DO 25 IDIM— 1 ,NDIM 
DO 25 ICE=1 ,NDIM 
DXDK(IDIM,ICE)=0. 

DO 24 KPE=1 , NPE 

DXDK(IDIM, ICE) = DXDK(IDIM, ICE) 

+ EX(KPE, IDIM)*DPHI(KPE, ICE) 

24 CONTINUE 

25 CONTINUE 
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GO TO (31,32,33), NDIM 

31 CONTINUE 
DETJB=DXDK(1 , 1) 

GO TO 37 

32 CONTINUE 

DETJ B=DXDK (1,1) *DXDK (2,2)- DXDK (1,2) *DXDK (2,1) 

GO TO 37 

33 CONTINUE 

DETJB = DXDK(1,1)*DXDK(2,2)*DXDK(3,3) 

+ DXDK(1 , 2)*DXDK(2 , 3)*DXDK(3 , 1) 

+ DXDK(2 , 1)*DXDK(3 , 2)*DXDK(1 , 3) 

- DXDK(1,1)*DXDK(2,3)*DXDK(3,2) 

- DXDK(2 , 2)*DXDK(1 , 3)*DXDK(3 , 1) 

- DXDK(3,3)*DXDK(2,1)*DXDK(1,2) 

C 

37 CONTINUE 

IF (DETJB. LE. 1 . E-15) THEN 

WRITE (6, 610) IELEM, IELF,NPE,NPRE 
DO 17 KPE=1 , NPE 

17 WRITE(6 , 615) KPE .NODEL(KPE) , (EX(KPE , IDIM) , IDIM=1 ,NDIM) 
WRITE(6 , 620) (PHI (K) , K=1 , NPE) 

DO 18 IDIM—1 , NDIM 

WRITE(6 , 625) (DPHI(K, IDIM) ,K=1,NPE) 

18 CONTINUE 

WRITE(6 ,630) DETJB ,( (DXDK(I ,J) ,J=1,NDIM) ,1=1, NDIM) 

STOP 
END IF 

610 FORMAT (2X, 'PROGRAM RUN TERMINATED AT SUB-ISOPEL DUE TO ', 
'SMALL DETJB' , /4X, ' IELEM=' , 15 , 2X, ' IELF-' ,12, 

2X, 'NPE=' ,12, 2X/NPRE-' ,12) 

615 FORMAT(3X, 'IPE ==' ,12, 2X, ' INODE-' , 15 , 2X, 'XDAT=' , 3E12 .4) 
620 FORMAT (4X, 'PHI =',5F10.4, 5(/5X,'PHI =',5F10.4)) 

625 FORMAT ( 4X, 'DPHI=' , 5F10 . 4, 5 (/4X, ' DPHI=' , 5F10 . 4) ) 

630 FORMAT ( 2X , ' DETJB=' , Ell . 4 , /2X , ' DXDK= ' , 3(/5X, 3E12 . 4) ) 

C 

GO TO (41,42,43), NDIM 

41 CONTINUE 
DKDX(1,1)=1. /DETJB 
GO TO 45 

42 CONTINUE 

DKDX(1 , 1)=DXDK(2 , 2) /DETJB 
DKDX(1,2)=- DXDK( 1,2) /DETJB 
DKDX(2 , 1) — DXDK(2 , 1) /DETJB 
DKDX (2,2)= DXDK (1,1) /DETJB 
GO TO 45 
C 

43 CONTINUE 
WRITE(6 , 635) NDIM 
STOP 

635 FORMAT (2X, 'TERMINATED IN SUB-ISOPEL FOR NDIM-', 13) 

C 

45 CONTINUE 
C 

C CALCULATE GLOBAL DERIVATIVES 



c 


DO 47 I D-IM-1 ,NDIM 

DO 47 IPE=T „NPE 

DPHX(IPE,IDIM)=0. 0 

DO 47 ICE=1,NDIM 

DPHX(IPE, IDIM) = DPHX(IPE.IDIM) 

+ DPHI (IPE , ICE) *DKDX(sICE , IDIM) 

47 CONTINUE 

GO TO ( 50',, 50*,, 5 © i; 5 0 , 5 2 , 53, 50>,.50 , 50 , 50',. 5O\,5O V 50'„5O>,.5O) 

I FLOW 

50' continue: 

WRITE (6, 6 50) I FLOW 
STOP 

6'50 ! FORMAT ( 2X "TERMINATED 1 AT SUB-ISOPEL FOR IFLOW=\I5) 

51 CALL. SHAE01S (PS ; I ,,DESI.„CXKS:,,NPRE)i 
GO TO 85 

52 CONTINUE 

DO 61 KDIM-1 , NDIM 
XGS(KDTM)-0. 

DO 61 KPE=1 , NPE 

XGS.(KDIM)=XGS (KDIM)+EX( KPE ,KDIM)*PHI(KPE) 

61 CONTINUE. 

ESI ('!.):= 1:.. 

DO 62 KDIM-1 ,.NDIM 

62 PSI(KDIM+l)=XGS(KDIM i )' 

GOi TO: 85 


53 CALL SHAP02 (PSI ,.DPSI , CXKS , NPRE) 

85 CONTINUE 

IF(IAXSY.EQ.l) THEN: 

RADUS=0 . 

DO 90 KEE=1 ,.NPE 
RADUS-RADUS+EX ( KPE,, 2)*PHI ( KPE) 
90 CONTINUE 
ENDIF 


RETURN’ 

END: 

c: 

Q^-M**.***ie^ieic^ie*.i(irk-k.Z*******^3*^***'^*U*'^****** : *5*********6^ 
SUBROUTINE LSHPI(ELK„DELK.,.S); 

G-X- IMPLICIT REAL*#- (A-H,0-Z> 

DIMENSION PLK(,3> ,.DPLK(.3> 

C: 

PLK(I) = (1 .. - S)/2 . 

PLK(.2 ) = (I . +S)/2 . 

DPLK(l) = -0 : . 5 
DPLK(2)= 0.5 
RETURN. 

END' 


4 . 4 : 



c 

C********l*********2*********3*********4*********5*********6 
SUBROUTINE LSHP2(PNK,DPNK,S) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

DIMENSION PNK(3) ,DPNK(3) 

C 

C 1-D QUADRATIC ELEMENT (NE= 123) 

C (S — 1. 0. 1.) 

PNK(1)=S*(S - 1 . )/2 . 

PNK(2)=1 . -S**2 
PNK(3)=S*(l.+S)/2. 

C 

DPNK( l)=S-0 . 5 
DPNK(2)=-2.*S 
DPNK(3)=S+0 . 5 
RETURN 
END 
C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE SHAP01 ( SHP , DSHP , CXKS , NPE) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION SHP(27) , DSHP (27 , 3) ,CXKS(3) 

C 

SHP(1)=1 . 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE SHAP02 ( SHP , DSHP , CXKS , NPE) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

DIMENSION SHP (27) , DSHP(27 , 3) , CXKS ( 3) 

SHP(1)=0. 333333333+0. 816496582*CXKS(2) 

SHP(2)=0. 333333333 -0 . 707106781*CXKS(1) -0 . 408248291*CXKS (2) 
SHP(3)=0 . 333333333+0 . 707106781*CXKS(1) -0 . 408248291*CXKS (2) 
RETURN 
END 
C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE SHAP23 (SHP , DSHP , CXKS , NPE) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION SHP(27) , DSHP (27 , 3) , PNK(3) , 

DPNK( 3) , PNE(3) , DPNE(3) , INDK(9) , INDE(9) , CXKS (3) 
DATA (INDK(KPE) ,KPE=1, 9) /1,2,3, 3,3,2, 1,1,2/, 

(INDE(KPE) ,KPE=1,9) /1,1,1, 2,3,3, 3,2,2/ 

C 7 6 5 

C 9 NODE QUADRATIC ELEMENT 894 

C 12 3 

CALL LSHP2(PNK, DPNK, CXKS (1) ) 

CALL LSHP2 ( PNE , DPNE , CXKS ( 2 ) ) 

C 

DO 10 KPE=1 ,NPE 
IPE=INDK(KPE) 

JPE=INDE(KPE) 

SHP(KPE)=PNK(IPE)*PNE(JPE) 
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DSHP(KPE , 1)=DPNK(IPE)*PNE(JPE) 

DSHP(KPE , 2)=PNK(IPE)*DPNE(JPE) 

10 CONTINUE 
RETURN 
END 
C 

C********l*********2*********3*********4*********5*********6 
SUBROUTINE PROCES (MAXNOD , MAXELM , MAXDOF , MXFRON) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE , NELEM , NPE , NPRE , NDIM , NEDOF , IFLOW , 
IAXSY, I ELF 

COMMON /CFRON/ MFRONF 

COMMON /CGRID/ X(4227 , 3) .NODES (27 , 1027) 

COMMON /CPROB/ IA(10),IPLOT 

COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 

C 

CALL PFRONT (NODES , NNODE , NELEM , NPE , MAXELM) 

CALL SHPLIB 
C 

CALL S FLOW (MAXNOD, MAXELM, MAXDOF, MXFRON) 

CALL PFLOW(MAXNOD, MAXELM, MAXDOF) 

IF(IPLOT.GE.l) CALL PLSDAT (MAXNOD, MAXELM, MAXDOF) 

C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6 
SUBROUTINE PFRONT (NODES .NNODE, NELEM, NPE, MAXELM) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION NODES (27, MAXELM) 

FIND LAST APPEARENCE OF EACH NODE AT FIRST ITERATION ONLY 

DO 30 INODE-1 , NNODE 
LASTE-0 

DO 20 KELEM-1, NELEM 
DO 10 IPE-l.NPE 
INODP-ABS (NODES ( IPE , KELEM) ) 

IF(INODP.NE. INODE) GO TO 10 
LASTE-KELEM 
LASTN-IPE 
GO TO 20 
10 CONTINUE 
20 CONTINUE 

NODES ( LASTN , LASTE ) = - INODE 
30 CONTINUE 
C 

RETURN 

END 

C 

C********i*********2*********3*********4*********5*********6*** 
SUBROUTINE SHPLIB 
C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE, NELEM, NPE, NPRE, NDIM, NEDOF, IFLOW, 
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IAXSY, IELF 

COMMON /CELEM/ EX(27 , 3) , EA(27 , 10) , ENUT(64) , ENUTl(B) , 

NODEL(27) , IELEM 

COMMON /CGAUS/ EXKS(3 , 64) , EW(64) ,MGAUS 

COMMON /CSHAP/ APHI(27 , 64) , APHX(27 , 3 , 64) , APSI (8 , 64) , AREA(64) , 
ARADUS(64) 

COMMON /ESHAP/ PHI (27) , DPHX(27 , 3) , PSI (27) , DXDK(3 , 3) , 

DKDX( 3 , 3 ) , CXKS ( 3 ) , DETJB , RADUS , IGAUS , JGAUS , LGAUS 
C 

REWIND 2 
C 

DO 100 IELEM=1 , NELEM 
CALL EXDAT 
C 

DO 50 LGAUS=1 , MGAUS 
DO 10 KDIM=1 , NDIM 
10 CXKS (KDIM)=EXKS(KDIM, LGAUS) 

CALL ISOPEL( IFLOW , IELF , IAXSY , NPE , NPRE , NDIM) 

DO 30 KPE=1 , NPE 

APHI (KPE , LGAUS )=PHI (KPE) 

DO 20 KDIM=1 , NDIM 

APHX(KPE , KDIM, LGAUS ) =DPHX(KPE , KDIM) 

20 CONTINUE 
30 CONTINUE 

IF(NPRE.GT.O) THEN 
DO 40 KPRE=1 , NPRE 
APSI (KPRE , LGAUS)=PSI (KPRE) 

. 40 CONTINUE 
END IF 

AREA ( LGAUS ) -»DETJ B*EW ( LGAUS ) 

IF ( IAXSY . EQ . 1 ) AREA ( LGAUS ) =RADUS*AREA ( LGAUS ) 

ARADUS ( LGAUS ) =RADUS 
50 CONTINUE 
C 

DO 90 L=l, MGAUS 

WRITE(2) (APHI (KPE, L) ,KPE=1,NPE) , ( ( APHX(KPE , K , L) ,KPE=1,NPE) , 
K=l, NDIM) .AREA (L) 

IF (NPRE . GT . 0) WRITE(2) (APSI (KPRE , L) ,KPkE=1 , NPRE) 

90 CONTINUE 

IF (IAXSY. EQ. 1) WRITE(2) (ARADUS (LGAUS ) , LGAUS=1 , MGAUS ) 

100 CONTINUE 
RETURN 
C 

C--- 1 2 3 4 5 6--- 

ENTRY SHPDAT 
DO 95 L=l, MGAUS 

READ ( 2 ) (APHI (KPE , L) ,KPE=1, NPE) , ( (APHX(KPE, K, L) ,KPE=1, NPE) , 
K-l.NDIM) ,AREA(L) 

IF (NPRE . GT . 0) READ ( 2 ) ( APSI (KPRE , L) ,KPRE=1, NPRE) 

95 CONTINUE 

IF(IAXSY.EQ.l) READ ( 2 ) (ARADUS (LGAUS) , LGAUS=1 , MGAUS) 

RETURN 

END 



C********l*********2*********3*********4*********5*********6*** 


SUBROUTINE S 1FLOW (NODES , NNODE , NELEM , NPE , 

NEDOF , NTDOF , I FLOW , MAXNOD , MAXELM , MAXDOF) 
C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDOF/ Al(11527) , IDBC(11527) , LDOF(4227) , L1D0F(4227) 
DIMENSION N0DES(27 .MAXELM) ,LIBDOF( 15) , LFLEL(27 , 15) 


C 

DATA (LIBDOF(IFL) , IFL=1 , 15) 

DATA (LFLEL(KPE , 5 ) , XPE=1 , 9 ) 
(LFLEL(KPE, 6) ,KPE=1,9) 
C 

DO 10 KELEM=1 , NELEM 
DO 10 KPE =1 ,NPE 
LDOF( ABS (NODES (KPE , KELEM) ) ) 
10 CONTINUE 
C 


/ 0, 0, 0, 0,18, 18, 0, 0, 0, 0, 

0 , 0 , 0 , 0 , 0 / 

/ 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 /, 

/ 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 / 


= LFLEL(KPE , I FLOW) 


LlDOF(l)=l 

DO 30 INODE==2 , NNODE 

LlDOF ( INODE ) =L1D0F ( INODE - 1 ) +LD0F ( INODE - 1 ) 

30 CONTINUE 

NEDOF=LIBDOF(IFLOW) 

NTDOF^LIDOF ( NNODE ) +LDOF ( NNOD E ) - 1 

IF(MAXDOF.GE. NTDOF) RETURN 

WRITE (6, 610) IFLOW.IPROB, MAXDOF, NTDOF 

STOP 

610 FORMAT (2X, 'TERMINATED IN SUB-S1FL0W FOR IFLOW=',I2, 

2X, 'IPROB=' ,12, 2X, 'MAXDOF=' ,16, 2X, ' NTDOF=' ,16) 

END 


C 

C********i*********2*********3*********4*********5*********6*** 
SUBROUTINE SEQVFL( IFLOW , MFRON , NTDOF , NDBC , NDIM , NNODE , 

MAXNOD , MAXELM , MAXDOF) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDOF/ Al(11527) , IDBC(11527) , LDOF(4227) , L1D0F(4227) 
COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 

COMMON /CFRON/ MFRONF 

COMMON /CINDX/ INDXF(27 , 3 , 15) , INDXP(27 , 15) 

COMMON /CPRS/ PELEM(4 ,1027) , PBCDAT , IPNOD(2) , IPDOF 
C 

MFRON=MFRONF 

C 

DO 10 KDOF=l, NTDOF 
IDBC(KDOF)=0 
Al(KDOF) =0. 

10 CONTINUE 
C 

DO 20 KNODE=l , NNODE 
DO 15 KDIM=1 ,NDIM 
IF(IBCA(KNODE,KDIM) .EQ.O) GO TO 15 
KD0F=L1D0F(KN0DE) - 1+KDIM 
IDBC(KDOF)= IBCA(KNODE.KDIM) 

Al(KDOF) = ADBC (KNODE , KDIM) 
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15 CONTINUE 
20 CONTINUE 
C 

KDBC=0 

DO 30 IDOF=l ,NTDOF 

IF(ABS(IDBC(IDOF)) .NE.O) KDBC=KDBC+1 
30 CONTINUE 
NDBC=KDBC 
C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6 
SUBROUTINE SFLOW(MAXNOD , MAXELM , MAXDOF , MXFRON ) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE , NELEM , NPE , NPRE , NDIM , NEDOF , I FLOW , 
IAXSY, I ELF 

COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 
COMMON /CFRON/ MFRONF 

COMMON /CGRID/ X(4227 , 3) , NODES (27 , 1027) 

COMMON /CITER/ CNVCF(IO) , ERROF(IO) ,RELAX(10) , ITERE, MAXIT 
COMMON /CPROB/ IA(10),IPL0T 
C 

ITERE=0 
1001 CONTINUE 

ITERE=ITERE+1 

I F ( ITERE . GT . MAXIT ) GO TO 2001 
C 

CALL FRONTS (NODES , IELIB , NNODE , NELEM , NEDOF , NPE , NPRE , 

NDIM , I FLOW , IAXSY , I ELF , MXFRON , MAXNOD , MAXELM , 
MAXDOF) 

C 

CALL SPRS4(A(1, 4) ,IBCA(1, 4) ,ERR0F(4) .MAXNOD, MAXELM, 
MAXDOF) 

C 

DO 120 IPROB-1,3 

IF(ERROF(IPROB) .GT.CNVCF(IPROB)) GO TO 1001 
120 CONTINUE 

IF(ERR0F(4) . GT.CNVCF(4) ) GO TO 1001 
C 

2001 CONTINUE 

CALL SPRS4(A(1 ,4) , IBCA(1 ,4) ,ERROF(4) , MAXNOD , MAXELM , 
MAXDOF) 

IF( ITERE. LE. MAXIT) RETURN 

CALL PLSDAT (MAXNOD, MAXELM, MAXDOF) 

C 

WRITE(6 , 688) MAXIT, ITERE 
688 FORMAT (2X, 'SOLUTION HAS FAILED TO CONVERGE', 

/4X, 'MAXIT=' ,15, 2X, 'ITERE=' ,15) 

STOP 

END 

C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE FRONTS (NODES , IELIB , NNODE , NELEM , NEDOF , NPE , NPRE , 



non o o 


NDIM, I FLOW , IAXSY , I ELF , MXFRON , MAXNOD , MAXELM , 
MAXDOF) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDOF/ Al(11527) , IDBC(11527) , LDOF(4227) , L1D0F(4227) 
COMMON /CELEM/ EX(27 , 3) , EA(27 , 10) , ENUT(64) , ENUT1 (8) 
NODEL(27) , IELEM 

DIMENSION NODES (27 .MAXELM) ,GK(167,167) ,GF(11527) , 

PNORM(167) ,LHEAD(167) ,EK(85,85) ,EF(85) , 

LOCEL(85) ,NDEST(85) 

C 

CALL SI FLOW ( NODES, NNODE , NELEM , NPE , 

NEDOF , NTDOF , IFLOW , MAXNOD , MAXELM , MAXDOF) 
CALL SEQVFL( IFLOW , MFRON , NTDOF , NDBC , NDIM , NNODE , MAXNOD , 

MAXELM .MAXDOF) 

C 

REWIND 1 
REWIND 2 
C 

C INITIALIZE HEADING AND GRAND FLUID MATRIX 

C 

NCRIT=MFRON- NEDOF 
NFRON=0 

DO 10 JFRON=l, MFRON 
DO 10 IFRON=l .MFRON 
GK ( I FRON , J FRON ) =0 . 

10 CONTINUE 

DO 20 IDOF=l, MAXDOF 
GF(IDOF)=0. 

20 CONTINUE 
C 

IELEM=0 
30 CONTINUE 

IELEM=IELEM+1 

CALL ELEMFL(EK . EF , NPE , NPRE , NDIM , NEDOF , IFLOW , 

IAXSY . IELF . MAXNOD , MAXELM , MAXDOF) 

C 

CREATE GLOBAL DOF ARRAY FOR EACH ELEMENT DOF 
IDOF=0 

DO 70 IPE-l.NPE 
INODE=NODES ( IPE , IELEM) 

N1D0F=L1D0F ( I ABS ( INODE ) ) 

NDOF=LDOF ( IABS ( INODE) ) 

DO 70 KDOF=l , NDOF 
IDOF=IDOF+l 

LOC EL ( I DOF ) =N 1 DOF+KDO F - 1 
IF ( INODE. LT.O) LOCEL(IDOF)=-LOCEL(IDOF) 

70 CONTINUE 

CONTRACT D.B.C. FOR ELEMENT SYSTEM OF EQUATIONS 

KDOF = 0 
NEWDOF= NEDOF 
DO 90 KDUM=1, NEDOF 
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KDOF = KDOF + 1 

IEQ = ABS (LOCEL(KDOF) ) 

IF(IDBC(IEQ) .EQ.O) GO TO 90 
C 

IF(KDOF.EQ.l) GO TO 81 
DO 80 ID0F=1 ,KDOF- 1 

EF(IDOF) = EF(IDOF) - EK(IDOF, KDOF) *Al( IEQ) 

IF (KDOF . EQ , NEWDOF) GO TO 80 
DO 71 JDOF=KDOF+l, NEWDOF 
EK(IDOF, JDOF-l)=EK(IDOF, JDOF) 

71 CONTINUE 

80 CONTINUE 
C 

81 CONTINUE 

IF(KDOF.EQ. NEWDOF) GO TO 86 
DO 85 IDOF-KDOF+1 , NEWDOF 

EF(IDOF-l) = EF(IDOF) -EK(ID0F,KD0F)*A1(IEQ) 

IF(KDOF.EQ.l) GO TO 83 

DO 82 JDOF=l , KDOF - 1 

EK( IDOF- 1 , JDOF) = EK(IDOF, JDOF) 

82 CONTINUE 
C 

83 CONTINUE 

DO 84 JDOF=KDOF+l , NEWDOF 

EK( IDOF - 1 , JDOF - 1 ) = EK(IDOF, JDOF) 

84 CONTINUE 

85 CONTINUE 
C 

86 CONTINUE 

DO 87 IDOF=l, NEWDOF 
EK( IDOF, NEWDOF) = 0. 

EK (NEWDOF , IDOF) - 0. 

EF (NEWDOF) =0. 

87 CONTINUE 

IF(KDOF.EQ. NEWDOF) GO TO 89 
DO 88 IDOFHKDOF+1, NEWDOF 
LOCEL(IDOF-l) = LOCEL(IDOF) 

88 CONTINUE 
C 

89 CONTINUE 

KDOF = KDOF - 1 
NEWDOF = NEWDOF - 1 

90 CONTINUE 
C 

C FIT EACH DOF INTO THE FRONT WIDTH EXTENDING IF NECESSARY 
C 

DO 120 IDOF=l, NEWDOF 
IEQ=LOCEL(IDOF) 

IF(NFRON.EQ.O) GO TO 95 
DO 94 IFRON=l , NFRON 
KFRON=IFRON 

IF(IABS(IEQ) . EQ . IABS (LHEAD(KFRON) ) ) GO TO 110 

94 CONTINUE 

95 CONTINUE 
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NFR0N=NFR0N+1 

IF(NFRON.LE.MFRON) GO TO 100 

WRITE (6,637) MXFRON , MFRON , NFRON , NCRIT , IELEM 

WRITE(6 ,638) ( LHEAD ( KFRON ) ,KFRON=l , NFRON) 

637 FORMAT (/2X, 'SUB -FRONTS --- FRONT WIDTH TO SMALL', 

- /4X, 'MXFRON-' ,15, 2X , ' MFRON= ',15, 2X, ' NFRON=' , 15 , 

- /4X, ' NCRIT- ' ,15, 2X, 'IELEM=' ,15, 

- /2X, ' LIST OF LHEAD DATA') 

638 FORMAT (2X, 1015) 

STOP 

100 CONTINUE 

NDEST(IDOF)=NFRON 
LHEAD (NFRON) =IEQ 
GO TO 120 
110 CONTINUE 

NDEST ( I DOF) =KFRON 
LHEAD (KFRON) =IEQ 
120 CONTINUE 


ASSEMBLE AN ELEMENT SYS. OF EQS . INTO A GLOBAL SYS. EQS 

DO 130 IDOF=l , NEWDOF 
IEQ=ABS (LOCEL(IDOF) ) 

GF(IEQ)=GF(IEQ)+EF(IDOF) 

IFRON=NDEST(IDOF) 

DO 130 JDOF=l, NEWDOF 
JFRON=NDEST ( JDOF) 

GK ( J FRON , I FRON ) =GK ( J FRON , IFRON)+EK(JDOF, IDOF) 

130 CONTINUE 

IF(NFRON. LT. NCRIT. AND. IELEM. LT.NELEM) GO TO 30 


CHECK THE LAST APPEARENCE OF EACH DOF 

140 CONTINUE 
PIV0T=0 . 0 

DO 170 IFRON=l , NFRON 

IF (LHEAD (I FRON) .GE.O) GO TO 170 

PIVOG=GK(IFRON, I FRON) 

IF(ABS(PIVOG) .LT.ABS (PIVOT)) GO TO 170 

PIVOT=PIVOG 

LPIVOT=IFRON 

170 CONTINUE 

IEQ=IABS (LHEAD (LPIVOT) ) 

IF (ABS (PIVOT) . GT . 1 . E- 10) GO TO 180 

WP.ITE (6,650) I PROB , I EQ , PIVOT , NCRIT , NFRON , I ELEM 
DO 171 IEQ-1, NEWDOF 
WRITE(6 , 652) IEQ.EF(IEQ) 

WRITE(6 , 654) (EK(IEQ.JEQ) ,JEQ-1 .NEWDOF) 

171 CONTINUE 
WRITE(6 ,656) 

WRITE(6 , 657) (NDEST ( JDOF) ,JDOF=l .NEWDOF) 
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WRITE(6 ,658) (LHEAD(IFRON) , IFRON=l , NFRON) 

WRITE (6, 659) LPIVOT , NFRON, GF(LPIVOT) 

WRITE(6 , 654) (GK(LPIVOT, JFRON) , JFRON=l , NFRON) 

WRITE(6 ,654) (GK(IFRON , LPIVOT) , IFRON=l .NFRON) 

STOP 

650 FORMAT (/2X, 'PROGRAM TERMINATED --- ILL-CONDITIONED MATRIX' , 
/4X, 'IPROB=' ,12, 2X, ' IEQ=' ,16, 2X, ' PIVOT=' , E12 . 4 , 

/4X, 'NCRIT=' ,15, 2X, 'NFRON=' ,15, 2X, ' IELEM=' , 15 , 

/4X, 'CURRENT ELEMENT IN PROCESS IELEM=',I5) 

652 F0RMAT(4X, 'IEQ=' ,12, 2X, ' EF(IEQ)=' , E12 . 4 , 2X, ' EK-DATA' ) 

654 F0RMAT(4X, 5E12 . 4) 

656 FORMAT (2X, 'CURRENT DATA IN THE GLOBAL MATRIX') 

657 FORMAT ( 2X , 'NDEST-DATA' , 20(/4X, 2013) ) 

658 FORMAT (2X, 'LHEAD- DATA' , 25 (/4X, 1016) ) 

659 FORMAT(2X, 'LPIVOT=' ,16, 2X, ' NFRON=' , 15 , 2X, ' GF=' , E12 . 4 , 

/2X, 'LIST OF PIVOTAL ROW AND COLUMN’) 

C 

180 CONTINUE 
C 

DO 190 IFRON=l, NFRON 
PNORM(IFRON)=GK(LPIVOT , IFRON) /PIVOT 
190 CONTINUE 

RHSID=GF( IEQ) /PIVOT 
GF(IEQ)=RHSID 
C 

IF(LPIVOT.EQ.l) GO TO 250 
DO 240 IFRON=l, LPIVOT- 1 
FACTOR=GK( IFRON, LPIVOT) 

C 

C UNDERFLOW MAY OCCUR IN THE FOLLOWING DO-200-LOOP IF 
C FACTOR IS SMALL. THE FOLLOWING STATEMENT NEED TO BE 
C CHANGED FOR DIFFERENT COMPUTERS. 

C 

DO 200 JFRON=l, LPIVOT -1 

GK ( I FRON , J FRON ) =GK ( I FRON , J FRON ) - FACTOR*PNORM ( J FRON ) 

200 CONTINUE 
C 

210 CONTINUE 

IF(LPIVOT.EQ. NFRON) GO TO 230 
DO 220 JFRON=LPIVOT+l, NFRON 

GK ( I FRON , J FRON - 1 ) =GK ( I FRON , J FRON ) - FACTOR* PNORM ( J FRON ) 

220 CONTINUE 
230 CONTINUE 

ITOTV=IABS( LHEAD ( IFRON)) 

GF ( ITOTV) =GF ( ITOTV) - FACTOR*RHS I D 
240 CONTINUE 
C 

250 CONTINUE 

I F ( LP I VOT . EQ . NFRON ) GO TO 300 
DO 290 IFRON=LPIVOT+l , NFRON 
FACTOR=GK( IFRON , LPIVOT) 

IF (LPIVOT. EQ. 1) GO TO 270 
DO 260 JFRON=l, LPIVOT -1 

GK( IFRON- 1 , JFRON )=GK (IFRON , JFRON) - FACTOR*PNORM( JFRON) 
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260 CONTINUE 
270 CONTINUE 

DO 280 JFRON=LPIVOT+1 , NFRON 

GK ( I FRON - 1 , J FRON - 1 ) =GK ( I FRON , J FRON ) - FACTOR*PNORM ( J FRON ) 

280 CONTINUE 

ITOTV=I ABS ( LHEAD ( I FRON ) ) 

GF ( ITOTV ) =GF ( ITOTV ) - FACTOR*RHS ID 
290 CONTINUE 
300 CONTINUE 

WRITE OUT NON -FIXED PIVOTAL EQUATION ON TAPE 

WRITE(l) NFRON, LPIVOT, ( LHEAD (I FRON ), PNORM( I FRON) , IFRON=l , NFRON) 

DO 320 IFR0N=1, NFRON 
GK ( I FRON , NFRON ) =0 . 0 
GK (NFRON, I FRON) =0 . 0 
320 CONTINUE 

IF (LPIVOT. EQ.^ NFRON) GO TO 340 
DO 330 IFRON=LPIVOT , NFRON - 1 
LHEAD (IFRON) =LHEAD( IFRON+1 ) 

330 CONTINUE 
340 CONTINUE 

NFRON=NFRON - 1 

ASSEMBLE, ELIMINATE, OR BACK- SUBSTITUTION 

IF(NFRON.GT.NCRIT) GO TO 140 
IF ( IELEM . LT. NELEM) GO TO 30 
IF (NFRON. GT.O) GO TO 140 

BACK- SUBSTITUTION 

DO 370 ITOTV=l , NTDOF -NDBC 
BACKSPACE 1 

READ(l) NFRON, LPIVOT, (LHEAD(IFRON) , PNORM(IFRON) , IFRON-1 , NFRON) 
IEQ=IABS (LHEAD (LPIVOT).) 

TEMPR=0 . 0 

PNORM ( LPIVOT )=0 . 0 

DO 360 IFRON=l, NFRON 

TEMPR=TEMPR - PNORM ( I FRON) *A1 ( I ABS ( LHEAD ( I FRON ) ) .) 

360 CONTINUE 

Al(IEQ) =GF( I EQ)+TEMPR 

BACKSPACE 1 
370 CONTINUE 

CALL SCNVFL (NODES , NNODE , NELEM, NPE , NPRE , NDIM , I FLOW , 

MAXNOD , MAXELM , MAXDOF) 

RETURN 

END 


C********l**^**^,*2******-***3*********4*.**-3Ht***3!r5***'******6*** 
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SUBROUTINE ELEMFL ( EK , EF , NPE , NPRE , NDIM , NEDOF , I FLOW , IAXSY , 

I ELF , MAXNOD , MAXELM , MAXDOF) 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CELEM/ EX(27 , 3) , EA(27 , 10) , ENUT(64) , ENUT1 (8) , 

NODEL(27) ,IELEM 

COMMON /CGAUS/ EXKS(3 , 64) , EW(64) ,MGAUS 
COMMON /CPNLT/ CXKPN(4 ,3,15), CWPN(4 , 15) , PTNUM 
COMMON /CINDX/ INDXF(27 , 3 , 15) , INDXP(27 , 15) 

COMMON /CMATE/ BFX(3) , DENSY.VISCY, PECLET 

COMMON /CSHAP/ APHI (27 , 64) , APHX(27 , 3 , 64) , APSI ( 8 , 64) ,AREA(64) , 
ARADUS(64) 

COMMON /CUSE2/ XK(3) ,XKN(3 ,3) ,XC(10) ,XB(3) ,XF(3) , PROD ,XWALL 
COMMON /CWIND/ WHI(27) ,DWHX(27 , 3) 

DIMENSION EK(85 ,85) , EF(85) , EQ(81 ,4) , QE(4 , 81) , EPM(4 , 4) , 

EPINV (4,4) , PIQE(4 , 81) , IROW(3) , JCOL(3) , DNUM(3) , 
DIFFU(3 , 3) 

C 

DO 1 IDOF-1, NEDOF 
EF(IDOF)=0 . 

DO 1 JDOF=l , NEDOF 
EK(IDOF, JDOF)=0 . 0 
1 CONTINUE 
C 

DO 3 KDOF-1 , NEDOF 
DO 3 KPRE=1 , NPRE 
EQ ( KDOF , KPRE ) =0 . 

QE ( KPRE , KDOF) =0 . 

3 CONTINUE 

DO 4 IPRE=1 , NPRE 
DO 4 JPRE=1 , NPRE 
EPM(IPRE , JPRE)=0 . 

4 CONTINUE 
C 

CALL ELMDAT 
CALL SHPDAT 
C 

DO 1000 LGAUS=1 , MGAUS 
XK(1)=VISCY 
DO 5 KDIM=1 , NDIM 
XC(KDIM)=0 . 

DO 5 KPE=1,NPE 

XC(KDIM)=XC(KDIM)+EA(KPE , KDIM)*APHI (KPE , LGAUS) 

5 CONTINUE 
C 

DO 7 KPE=1 ,NPE 

WHI (KPE)=APHI (KPE , LGAUS) 

7 CONTINUE 

DO 10 KDIM=1 , NDIM 
DO 10 KPE=1 ,NPE 

DWHX(KPE , KDIM)=APHX (KPE ,KDIM , LGAUS) 

10 CONTINUE 
C 

DO 30 IPE=1 ,NPE 
DO 11 KDIM=1 ,NDIM 
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IROW(KDIM)=INDXF(IPE,KDIM, IFLOW) 

EF( IROW(KDIM) )=EF( IROW(KDIM) )+WHI ( IPE) *BFX(KDIM) *AREA( LGAUS ) 

11 CONTINUE 

DO 25 JPE=1,NPE 
DO 12 KDIM=1,NDIM 
JCOL(KDIM)=INDXF(JPE,KDIM, IFLOW) 

12 CONTINUE 

CONVC=0 . 

DIFF=0 . 

DO 15 IDIM=1 ,NDIM 

CONVC=CONVC+WHI (IPE)*DENSY*XC ( IDIM)*APHX( JPE , IDIM , LGAUS) 
*AREA( LGAUS) 

DIFF HDIFF +DWHX(IPE,IDIM)*XK(1)*APHX(JPE, IDIM, LGAUS) 
*AREA( LGAUS) 

DO 14 JDIM=1 ,NDIM 

DIFFU( IDIM, JDIM)=DWHX(IPE,JDIM)*XK(1)*APHX( JPE, IDIM, LGAUS) 
*AREA( LGAUS) 

14 CONTINUE 

15 CONTINUE 

DO 20 IDIM=1 ,NDIM 

EK(IROW(IDIM) , JCOL(IDIM))=EK(IROW(IDIM) , JCOL(IDIM) )+CONVC 

+DIFF 

DO 19 JDIM=1 , NDIM 

EK(IROW(IDIM) , JCOL(JDIM) )=EK(IROW(IDlM) .JCOL(JDIM)) 

+DIFFU(IDIM , JDIM) 

19 CONTINUE 

20 CONTINUE 

IF(IAXSY.EQ.l) THEN 

THETV = 2 ,*XK(1)*WHI(IPE)*APHI(JPE, LGAUS) 

/ARADUS ( LGAUS ) **2*AREA ( LGAUS ) 

EK(IROW(2) ,JCOL(2))=EK(IROW(2) , JC0L(2) )+THETV 
END IF 

25 CONTINUE 
30 CONTINUE 

DO 52 IPRE=1,NPRE 
DO 52 JPRE=1,NPRE 

EPM ( I PRE , J PRE ) =EPM ( I PRE , J PRE ) +APSI ( I PRE , LGAU S ) 

*APSI (JPRE , LGAUS)*AREA(LGAUS) . 

52 CONTINUE 

DO 53 IPE=1 ,NPE 

IROW( 1)=INDXF( IPE, 1 , IFLOW) 

IROW ( 2 ) =INDXF ( I PE , 2 , IFLOW) 

DO 53 JPRE=1 ,NPRE 

EQ ( IROW ( 1 ) , JPRE)=EQ(IROW(l) , JPRE)+DWHX(IPE, 1) 

*APS I (J PRE , LGAUS ) *AREA (LGAUS ) 

EQ(IROW(2) , J PRE ) =EQ ( IROW ( 2 ) , JPRE)+DWHX(IPE, 2) 

*APS I (JPRE , LGAUS ) *AREA ( LGAUS ) 
IF(IAXSY.EQ.l) EQ ( IROW ( 2 ) , JPRE)=EQ(IROW(2) , JPRE)+WHI (IPE) 

*APS I (JPRE , LGAUS ) /ARADUS ( LGAUS ) *AREA ( LGAUS ) 

53 CONTINUE 
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DO 54 IPRE=1,NPRE 

DO 54 JPE=1 ,NPE 

J COL( 1 ) =INDXF( JPE , 1 , 1 FLOW) 

J COL ( 2 ) =INDXF ( J PE , 2 , I FLOW) 

QE(IPRE, JCOL(l) )=QE(IPRE , JCOL(l) ) +PTNUM*APS I ( I PRE , LGAU S ) 

*APHX ( J PE , 1 , LGAUS ) *AREA ( LGAUS ) 

QE ( I PRE , J COL ( 2 ) ) =QE ( I PRE , J COL ( 2 ) ) +PTNUM*AP S I ( I PRE , LGAU S ) 

*APHX(JPE, 2 , LGAUS )*AREA( LGAUS) 
IF(IAXSY.EQ.l) EK(IPRE, JCOL(2) )=EK(IPRE, JCOL(2) ) + PTNUM 

*APS I ( I PRE , LGAUS ) *APHI ( JPE , LGAUS ) *AREA ( LGAUS ) /ARADUS ( LGAUS ) 
54 CONTINUE 
1000 CONTINUE 
C 

DETM - EPM(1,1)*EPM(2,2)*EPM(3,3) 

+ EPM(1,2)*EPM(2,3)*EPM(3,1) 

+ EPM(2,1)*EPM(3,2)*EPM(1 1 3) 

- EPM(1,1)*EPM(2,3)*EPM(3,2) 

- EPM(2 , 2)*EPM(1 , 3)*EPM(3 , 1) 

- EPM(3,3)*EPM(2,1)*EPM(1,2) 

EPINV(1 , 1)=(EPM(2 , 2)*EPM(3 , 3) -EPM(3 , 2)*EPM(2 , 3) )/DETM 
EPINV ( 1 , 2 ) = ( EPM (1,3) *EPM (3,2) - EPM (1,2) *EPM (3,3)) /DETM 
EPINV(1 , 3)=(EPM(1 , 2)*EPM(2 , 3) -EPM(2 , 2)*EPM(1 , 3) ) /DETM 
EPINV (2 , 1)=(EPM(2 , 3)*EPM(3 , 1) -EPM (2 , 1)*EPM(3 , 3) )/DETM 
EPINV (2 , 2)=(EPM(1 , 1)*EPM(3 , 3) -EPM(3 , 1)*EPM(1 , 3) ) /DETM 
EPINV(2,3)=(EPM(2,1)*EPM(1,3)-EPM(1,1)*EPM(2,3))/DETM 
EPINV (3 , 1)=(EPM(2 , 1)*EPM(3 , 2) -EPM (2 , 2)*EPM(3 , 1) ) /DETM 
EPINV (3 , 2)=(EPM(1 , 2)*EPM(3 , 1) -EPM(1 , 1)*EPM(3 , 2) )/DETM 
EPINV(3,3)-(EPM(1,1)*EPM(2,2) -EPM(2 , 1)*EPM(1 , 2) )/DETM 
C 

DO 83 IPRE=1 , NPRE 
DO 83 JDOF-l.NEDOF 
PIQE(IPRE, JDOF)=0 . 

DO 82 KDUM—1 , NPRE 

PIQE ( IPRE , JDOF) =PIQE ( I PRE , JDOF) +EPINV ( I PRE , KDUM) 

*QE(KDUM, JDOF) 

82 CONTINUE 

83 CONTINUE 

DO 85 IDOF-l.NEDOF 
DO 85 J DOF-1, NEDOF 
DO 84 KDUM-1 , NPRE 

EK( IDOF , JDOF) =EK( IDOF , JDOF) +EQ ( IDOF , KDUM) *PIQE (KDUM , JDOF) 

84 CONTINUE 

85 CONTINUE 
C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6*** 
SUBROUTINE S CNVFL( NODES, NNODE , NELEM , NPE , NPRE, NDIM, I FLOW, 

MAXNOD , MAXELM , MAXDOF) 

C-X- IMPLICIT REAL*8 (A-H.O-Z) 

COMMON /CDOF/ Al(11527) , IDBC(11527) ,LDOF(4227) ,L1D0F(4227) 
COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 

COMMON /CITER/ CNVCF(IO) , ERROF(IO) , RELAX (10) , ITERE ,MAXIT 
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op o on Poo 


COMMON /CPRS/ PELEM(4 , 1027) , PBCDAT, IPN0D(2) , IPDOF 
DIMENSION NODES (27 .MAXELM) , KERR (4) ,DELA(4) 

Al NEW SOLUTION OBTAINED IN SUB-FRONTS 

DO l'K-1,4 

1 ERR0F(K)=0. 

AVELY=0 . 

DO 5 KNODE=l , NNODE 
IDOF=LlDOF (KNODE ) - 1 
ADUM=0. 

DO 2 KDIM=1 ,NDIM 
ADUM=ADUM+A1 ( IDOF+KDIM) **2 

2 CONTINUE 
ADUM=ADUM**0 . 5 

IF ( ADUM . GT . AVELY) AVELY=ADUM 
5 CONTINUE 

DO 10 KN0DE=1, NNODE 
ID0F=L1D0F (KNODE) -1 

DO 7 KDIM=1 ,NDIM 
KDOF=IDOF+KDIM 

DELA(KDIM)=ABS(A1(KD0F) -A(KNODE.KDIM) )/AVELY 
IF(DELA(KDIM) . GT. ERROF(KDIM) ) THEN 
ERROF(KDIM) =DELA (KDIM) 

KERR ( KD IM ) =KNODE 
ENDIF 
7 CONTINUE 
10 CONTINUE 

DO 15 KNODE=l, NNODE 
ID0F=L1D0F (KNODE) -1 
DO 12 KDIM=1 , NDIM 
KDOF=IDOF+KDIM 
IF(IDBC(KDOF) . EQ.l) THEN 
A(KN0DE , KDIM)=A1 (KDOF) 

ELSE 

A(KNODE ,KDIM) = (1 . -RELAX(KDIM) )*A(KNODE ,KDIM) 

+RELAX(KDIM) *Al (KDOF) 

ENDIF 

12 CONTINUE 
15 CONTINUE 


WRITE(6 , 630) ITERE, (K,KERR(K) , ERROF(K) ,K=1 ,NDIM) 

630 FORMAT (2X, 'SUB- SCNVFL ITERE=' ,15, 

4 (/4X , ' KDIM= ',12, 2X, 'N0DE=' ,15, 2X, ' ERROF=' , E12 . 4) ) 

C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6*** 



SUBROUTINE SPRS4 ( P , IBCP , PERROR , MAXNOD , MAXELM , MAXDOF) 

X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE , NELEM, NPE , NPRE , NDIM, NEDOF, I FLOW , 

IAXSY , I ELF 

COMMON /CELEM/ EX(27 , 3) , EA(27 , 10) ., ENUT(64) , ENUT1(8) , 
NODEL(27) , IELEM 

COMMON /CGAUS/ EXKS (3 , 64) , EW(64) .MGAUS 
COMMON /CLSCF/ XKSNOD(27 , 3 , 11) , TM(4 , 4) 

COMMON /CMATE/ BFX(3) , DENSY, VISCY, PECLET 
COMMON /CPNLT/ CXKPN(4 ,3,15), CWPN(4 , 15) , PTNUM 
COMMON /CPROB/ IA(10),IPLOT 

COMMON /CPRS/ PELEM(4 , 1027) , PBCDAT , IPNOD(2) , IPDOF 
COMMON /CSHAP/ APHI (27 , 64) , APHX(27 , 3 , 64) , APSI (8 , 64) , AREA(64) 
ARADUS(64) 

DIMENSION P (MAXNOD) , IBCP (MAXNOD) , EK(85 , 85) ,EF(85) ,EKI(4,4) , 
PNLT(4) ,DADX(3) ,DIV(3) ,CXKS(3) ,PSIZ(27) , 

DPSIZ(27 , 3) , PSIN0D(4 , 27) 


REWIND 2 
IPROF=0 

DO 8 KPE=1,NPE 

DO 1 KDIM=1 ,NDIM 

CXKS (KDIM)=XKSNOD (KPE , KDIM, IELF) 

1 CONTINUE 

GO TO (2,2,2,2,10, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2), IFLOW 

2 CONTINUE 
WRITE(6 , 605) IFLOW 
STOP 

605 FORMAT (2X, 'TERMINATED AT SUB- SPRS4 FOR IFLOW=',I5) 

3 CONTINUE 

CALL SHAP02(PSIZ,DPSIZ, CXKS, NDIM) 

DO 7 KPRE=1 , NPRE 

PS INOD ( KPRE , KPE ) =PS IZ ( KPRE ) 

7 CONTINUE 

8 CONTINUE 

10 CONTINUE 

DO 11 KNODE-1 , NNODE 
P(KNODE) =0. 

IBCP(KNODE)=0 

11 CONTINUE 

PERROR=0 . 

PMAX =0 . 

DO 13 KELEM=1, NELEM 
DO 12 KPRE=1 ,NPRE 
DUM=ABS ( PELEM ( KPRE , KELEM) ) 

IF(DUM.GT.PMAX) PMAX=DUM 

12 CONTINUE 

13 CONTINUE 


DO 1000 IELEM=1, NELEM 



DO 14 IPE=1 ,NPE 
EF(IPE)=0. 

DO 14 JPE=1,NPE 
EK(IPE, JPE)=0. 

14 CONTINUE 
CALL EXDAT 
CALL ELMDAT 
CALL SHPCAT 

DO 100 LGAUS=1 , MGAUS 
XF=0 . 

DO 16 KDIM=1,NDIM 
DADX(KDIM)=0. 

DO 15 KPE=1 , NPE 

DADX(KDIM)=DADX(KDIM)+EA(KPE , KDIM) *APHX (KPE ,KDIM , LGAUS ) 

15 CONTINUE 
XF=XF+D ADX ( KD IM ) 

16 CONTINUE 

IF(IAXSY.EQ.l) THEN 
WELY=0. 

DO 20 KPE=1 ,NPE 

WELY=WELY+EA ( KPE , 2 ) *APHI (KPE , LGAUS ) 

20 CONTINUE 

XF = XF + WELY/ARADUS ( LGAUS) 

END IF 

DO 25 IPRE=1 , NPRE 

EF ( I PRE ) =EF ( I PRE ) +APS I ( IPRE , LGAUS ) *PTNUM*XF*AREA ( LGAUS ) 

DO 25 JPRE=1 , NPRE 

EK ( I PRE , J PRE )=EK(I PRE , J PRE ) +APS I ( IPRE , LGAUS ) 

*APS I ( J PRE , LGAUS ) *AREA ( LGAUS ) 

25 CONTINUE 
100 CONTINUE 

DETM=EK(1 , 1)*EK(2 , 2)*EK(3 , 3) + EK(1 , 2)*EK(2 , 3)*EK(3 , 1) 
+EK(2,1)*EK(3,2)*EK(1 J 3) - EK(1, 1)*EK(2, 3)*EK(3 , 2) 
-EK(2,2)*EK(1,3)*EK(3,1) - EK(3 , 3)*EK(2 , 1)*EK(1 , 2) 

EKI (1 , 1)=(EK(2 , 2)*EK(3 , 3) -EK(3 , 2)*EK(2 , 3) )/DETM 
EKI (1 , 2)=(EK(1 , 3)*EK(3 , 2) -EK(1 , 2)*EK( 3 , 3) )/DETM 
EKI (1 , 3)=(EK(1 , 2)*EK( 2 , 3) -EK(2 , 2)*EK(1 , 3) )/DETM 
EKI (2 , 1)=(EK(2 , 3)*EK(3 , 1) -EK(2 , 1)*EK(3 , 3) )/DETM 
EKI (2 , 2)=(EK(1 , 1)*EK(3 , 3) -EK(3 , 1)*EK(1 , 3) )/DETM 
EKI (2 , 3)=(EK(2 , 1)*EK( 1 , 3) -EK(1 , 1)*EK(2 , 3) )/DETM 
EKI (3 , 1)=(EK(2 , 1)*EK(3 , 2) -EK(2 , 2)*EK(3 , 1) )/DETM 
EKI(3,2)=(EK(1,2)*EK(3,1)-EK(1,1)*EK(3,2))/DETM , 

EKI (3 , 3) = (EK(1 , 1)*EK(2 , 2) -EK(2 p 1)*EK(1 , 2) )/DETM 

DO 50 IPRE=1 , NPRE 
PNLT(IPRE)=0 . 

DO 40 JPRE=1 ,NPRE 

PNLT ( I PRE ) =PNLT ( I PRE ) - EKI(IPRE, JPRE)*EF(JPRE) 

40 CONTINUE 
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on on 


50 CONTINUE 
C 

IF(PMAX.LE.l.E-6) GO TO 61 
DO 60 KPRE=1 ,NPRE 

DUM=ABS (PNLT(KPRE) -PELEM(KPRE , IELEM) ) /PMAX 
IF (DUM . GT . PERROR) THEN 
PERROR=DUM 
KPRNOD=NODEL (NPE) 

ENDIF 

60 CONTINUE 
C 

61 CONTINUE 
DO 65 KPRE=1 , NPRE 
PELEM ( KPRE , I ELEM ) =PNLT (KPRE) 

65 CONTINUE 

GO TO (71,71,71,71,75, 76,71,71,71,71, 71,71,71,71,71) 

I FLOW 

71 CONTINUE 

WRITE(6 , 672) IFLOW 
STOP 

672 FORMAT (2X, 'TERMINATED AT SUB- SPRS4 FOR IFLOW=',I5) 

PRESSURE INTERPOLATION POLYNOMIALS OF THE FORM [1,X,Y] - 

75 CONTINUE 
DO 140 KPE=1 ,NPE 
KNODE=NODEL(KPE) 

IBCP(KNODE)-IBCP(KNODE)+l 
PDUM - PNLT(l) 

DO 135 KDIM=1 ,NDIM 

PDUM - PDUM + PNLT(KDIM+1)*EX(KPE,KDIM) 

135 CONTINUE 

P(KNODE) = P(KNODE) + PDUM 
140 CONTINUE 
GO TO 145 

--- NEW PRESSURE INTERPOLATION POLYNOMIALS - 

76 CONTINUE 
DO 144 KPE=1 , NPE 
DO 142 KDIM=1 , NDIM 
CXKS (KDIM)=XKSNOD(KPE , KDIM , I ELF) 

142 CONTINUE 
KNODE=NODEL(KPE) 

IBCP(KNODE)=IBCP(KNODE)+l 
PDUM=0 . 

DO 143 KPRE=1,NPRE 

PDUM = PDUM + PNLT (KPRE) *PSIN0D (KPRE, KPE) 

143 CONTINUE 
P(KNODE)=P(KNODE) + PDUM 

144 CONTINUE 

145 CONTINUE 
C 

1000 CONTINUE 

687 FORMAT (2X, 'SUB- SPRS IELEM=',I5, 2X, 'KPE=' , 12 , 



no. on non 


2X, 'KNODE=' ,16, /2X,'X=' ,E12.4, 2X, ' Y-' , E12 . 4 , 

2X, 'P=' ,E12.4) 

C 

WRITE (6, 6 50) KPRNOD , PERROR 

650 FORMAT (2X, 'KDIM=4 KPRN0D=',I6, 2X, ' PERR0R=' , E12 . 4) 

C 

DO 150 KN0DE=1 , NNODE 
P(KN0DE)=P(KN0DE) /FLOAT (IBCP(KNODE) ) 

150 CONTINUE 
C 

PREF=P ( I PNOD ( 2 ) ) 

DO 160 KNODE=l , NNODE 
P ( KNODE )=P( KNODE) -PREF+PBCDAT 
160 CONTINUE 
C 

RETURN 

END 


********i*********2*********3*********4*********5*********6*** 
SUBROUTINE PFLOW(MAXNOD , MAXELM , MAXDOF) 

-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE , NELEM ,NPE , NPRE , NDIM , NEDOF, IFLOW , 

IAXSY , I ELF 

COMMON /CGRID/ X(4227 , 3) ,NODES(27 , 1027) 

COMMON /CITER/ CNVCF(IO) , ERROF(IO) .RELAX (10) , ITERE ,MAXIT 
COMMON /CPROB/ IA(10),IPLOT 

COMMON /CPRS/ PELEM(4 , 1027) , PBCDAT, IPNOD(2) , IPDOF 
COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 

WRITE(6 , 650) ITERE, IFLOW 

650. FORMAT ( 2X , 'ENTRY-PFLOW ITERE=' ,14, 2X, 'IFLOW-' ,12) 

DO 50 KNODE-1, NNODE 

WRITE(6 , 660) KNODE, (A(KN0DE, IDUM) , IDUM=1 ,4) 

50 CONTINUE 
RETURN 

660 FORMAT (2X, 15 , 5E12 . 4) 

--1 2- 3 4 5 6--- 

ENTRY PLSDAT(MAXNOD, MAXELM, MAXDOF) 

WRITE(7 , 605) 

DO 10 KNODE-1, NNODE 

WRITE(7 , 606) KNODE, (X(KNODE.KDUM) ,KDUM=1 , 3) 

10 CONTINUE 

605 FORMAT ( 2X , ' KNODE X Y Z') 

606 FORMAT(2X, 15 , 2X, 3E16 . 8) 

WRITE(7 ,612) 

612 FORMAT (2X, 'NODE CONNECTIVITY DATA') 

DO 15 KELEM-1, NELEM 

WRITE(7 , 615) KELEM, (NODES (KPE.KELEM) ,KPE=1,NPE) 

15 CONTINUE 

615 FORMAT (2X, 15 , 2X, 1017 , 2(/9X,10I7)) 
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WRITE(7 , 620) 

DO 20 KN0DE=1 , NNODE 

WRITE(7 , 621) KNODE, (I BCA( KNODE, KPROB) ,KPR0B=1,3) 

20 CONTINUE 
C 

WRITE(7 , 624) 

WRITE(7 , 625) PBCDAT , (IPNOD(K) ,K=1 , 2) 

620 FORMAT ( 2X , 'IBC- DATA FOR IA=1 , 2 , 3 , 4 , 6 ' ) 

621 FORMAT (4X, 15 , 2X, 2013) 

624 FORMAT (2X, 'PBCDAT AND IPNOD(l-2)') 

625 FORMAT ( 2X, E14 . 6 , 2X,2I7) 

C 

WRITE(7 , 630) 

DO 30 KNODE=l, NNODE 

WRITE(7 , 631) KNODE, (A(KNODE.K) ,K=1,4) 

30 CONTINUE 

630 FORMAT (2X, ' KNODE U V W' , 

P' ) 

631 FORMAT (2X, 15 ,4E17 . 9) 

C 

WRITE(7 , 640) 

DO 40 KNODE=l, NNODE 

WRITE(7 , 641) KNODE, (ADBC (KNODE ,K) ,K=1,3) , ADBC(KNODE, 6) 

40 CONTINUE 

640 FORMAT (2X, 'ADBC -DATA FOR U, V, AND W' ) 

641 FORMAT(2X, 15 ,4E17 . 9) 

C 

RETURN 

END 

C 

C********l*********2*********3*********4*********5*********6 

SUBROUTINE USER 

C-X- IMPLICIT REAL* 8 (A-H.O-Z) 

COMMON /CDESC/ NNODE, NELEM,NPE,NPRE,NDIM,NEDOF, I FLOW, 

IAXSY, I ELF 

COMMON /CELEM/ EX(27 , 3) , EA(27 , 10) , ENUT(64) , ENUT1(8) , 

NODEL(27) ,IELEM 

COMMON /CGAUL/ CLXKS (4,4), CLW (4,4), NGAUS 
COMMON /CGAUS/ EXKS(3 , 64) , EW(64) ,MGAUS 
COMMON /CGRID/ X(4227 , 3) ,NODES(27 , 1027) 

COMMON /CMATE/ BFX(3) , DENSY, VISCY, PECLET 

COMMON /CSHAP/ APHI (27 , 64) , APHX(27 , 3 , 64) , APSI (8 , 64) , AREA(64) , 
ARADUS (64) 

COMMON /CUSE2/ XK(3) ,XKN(3 , 3) ,XC(10) ,XB(3) ,XF(3) , PROD ,XWALL 
COMMON /CFLOW/ A(4227 , 10) , ADBC(4227 , 10) , IBCA(4227 , 10) 
DIMENSION PNK(3) ,DPNK(3) 

C 

C 1 2 3 4 5- 6- - - 

ENTRY EXDAT 

DO 4 KPE=1 , NPE 

KNODE=ABS (NODES (KPE , IELEM) ) 

NODEL(KPE)=KNODE 
DO 3 KDIM=1 ,NDIM 
EX(KPE , KDIM)=X(KNODE ,KDIM) 
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3 CONTINUE 

4 CONTINUE 
RETURN 

C 

C 1- 2 3 4-- 5 6 

ENTRY ELMDAT 

DO 6 KPE-l.NPE 

KNODE=ABS (NODES (KPE , IELEM) ) 

NODEL(KPE)=KNODE 

DO 5 KPROB=l , 10 

EA ( KPE , KPROB ) =A ( KNODE , KPROB ) 

5 CONTINUE 

6 CONTINUE 
C 

RETURN 

END 

******************************** BOTTOM OF DATA ******************************* 
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APPENDIX II 


INPUT DATA FOR NSFLOW/P 


The required input data to solve the incompressible, laminar flows is described 
below. The computational sequence is controlled by the macro-instruction data [27] 
in the main program. These macro-instruction data are "INIT," "PREP," "PROC," 
"CONT," AND "END;" and these data have to start from the fifth column of each 
card. The function of these data are described below. 

"INIT" — Initialize dimensioned variables. 

"PREP" — Call the SUBROUTINE PREP to read in the descriptive data for each 
flow problem. 

"PROC" — Call the SUBROUTINE PROCES to solve the Navier-Stokes equations. 

"CONT" — Continue computation for the next flow problem. 

"END" — Terminate the computation. 

The descriptive data for a specific flow case are read into the computer program 
in the SUBROUTINE PREP. The sequence to read in the various descriptive data is 
also controlled by the macro-instruction data. The macro-instruction data used in the 
SUBROUTINE PREP are listed below. The function for each of these macro-instruction 
data and a set of specific data followed by each of these macro-instruction data are 
described below. The macro-instruction data used in the SUBROUTINE PREP have to 
start from the first column of each card. In most of the cases, a comment card has 
been used to clarify the input data to be prepared. 

1. "DESC" — Read in the general descriptive data. 

IFLOW — =5, Solve two-dimensional flows using the pressure interpola- 
tion polynomials of the form (l,x,y); = 6, Solve two-dimensional 
flows using the new pressure interpolation method. 

NDIM — Dimension of the problem. 

NGAUS — Number of Gauss points in each coordinate direction. (Ngaus 
= 3 has been tested). 

MFRONF — Frontal width. 

2. "CNTL" — Control parameters. 

NNODE — Number of nodes. 

. NELEM — Number of elements. 

IAXSY — = 0 for two-dimensional case, and = 1 for axisymmetric case. 

IPLOT — = 1 to write the computational results on a disk file. 


65 



3. "ELEM" — Call the SUBROUTINE ELEM to generate the node connectivity 

data. The input data for the subroutine is described below. 
Again, some of the data are followed by a comment card. 

NBLOC — Number of blocks to generate the node connectivity data. 

IEL1 — The first element number in each block. 

(NODES(IPE,IELl) ,IPE=1,NPE) — Node connectivity data for the first 

element in each block. NPE is the number of nodes in an element 

NEL(KDIM) — Number of elements in each coordinate direction. 

INCREL(KDIM) — Incremental element number in each coordinate direction. 

(INCNOD(K ,KDIM) ,K=1,NPE) — Increment of the connectivity data for 
each coordinate direction. 

4. "NODE" — Call the SUBROUTINE RNODE to generate the grid coordinate 

data. 

NBLOC — Number of blocks for the coordinate data generation. . 

METHOD — = 1, To read in the coordinate data on the physical domain; 

= 2, To read in the coordinate data on the computational 
element. In this case, isoparametric mapping is used for 
grid generation. 

Description of the input data for METHOD=l 


NODG1 — The first node number in each block. 

INCRX, INCRY, INCRZ — Incremental node numbers in each coordinate 
direction . 

NDAT — Number of grid points in each coordinate direction. 

(DELXQKE ,KDIM) ,IKE=1,NDAT) — An array of physical coordinate data 
in each coordinate direction. 

Description of input data for METHOD=2 


((XNOD(KPE,KDIM) ,KDIM=1,KPE=1,NPE) — Coordinate data of the block. 
The sequence of node numbers should be the same as that of 
the computational element. 

NODG1, INCRX, INCRY, INCRZ - The same as above. 

NDAT — The same as above. 

(DELX(IKE ,KDIM) ,IKE=1,NDAT) — An array of cordinate data defined on 
the computational element in each coordinate direction. 
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5. "MATE" — Material property data. 

VISCY — Molecular viscosity of the fluid. 

DENSY — Density of the fluid. 

(BFX(K) ,K=1,NDIM) — Body force in each coordinate direction. 

6. "ITER" — Iteration parameters. 

MAXIT — Maximum number of iterations. 

(RELAX(K) ,K=1, 10) — tinder -relaxation numbers; K = 1, 2, and 3 for 
the x- , y- , z- momentum equations, respectively; K = 4 for 
pressure; rest of these under- relaxation numbers are not used 
as yet . 

(CNVCF(K) ,K=1, 10) — Convergence criteria, use is the same as above. 

7. "IA##" - Call the SUBROUTINES RINIT and RBC1 to read in the initial 

guess and the boundary condition data for flow equations. 

(## = 1, 2, and 3 for u, v, and w, respectively; IA05 through 
I A 10 are not used as yet.) 

Input data for SUBROUTINES RINIT and RBC1 

NREC — Number of records. 

N1 — The first node number. 

N2 — The last node number. 

INCNOD — Incremental node number. 

ADATA — Real variable to be assigned as the initial guess. 

8. "IA04" — Input data for pressure. 

PBCDAT — A real variable for pressure boundary condition. 

IPNOD(l) — A pressure node number for which the pressure boundary 
condition is specified. 

IPNOD(2) — A velocity node number to prescribe a reference pressure. 

9. "INCL" — Include re-start data. 

10. "END" — Return the control to the main program. 
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A. 2.1 Cavity Flow for Re = 10,000 


********************************* TOP OF DATA *************************** 
CAVITY FLOW FOR REYNOLDS NUMBER=10000 (CAVT91) 

****INITIALIZE DIMENSIONED VARIABLES**** 

****PREPARE INPUT DATA **** 

DESCRIPTIVE DATA - - 

I FLOW, NDIM, NGAUS , MFRONF , 

6, 2, 3, 165, 

CNTL PARAMETERS - - 

NNODE , NELEM, IAXSY, IPLOT, 

4225, 1024, 0, 1, 

MATERIAL PROPERTY OF FLUID - 

VISCY, DENSY, BFX(l-2), 

0.0001225, 1.225, 0. , 0. , 

ITERATION PARAMETERS - - 

MAXIT, RELAX(l-lO) , CNVCF(l-lO) 


100 , 

0.8, 0. 

8, 1., 

1 . , 

1., - 


1 ., 1 . 

1. , 

1- , 

1-, 


l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 

l.E-4, 


NODE COORDINATE DATA - GRID GENERATION 
NUMBER OF BLOCKS (NBLOC) 

1 , 

GRID GENERATION METHOD FOR IBLOC=l (METHOD) 

2 , 

NODE COORDINATE DATA ( (XNOD(KPE.KDIM) ,KDIM=1 ,NDIM) ,KPE=1 ,NPE) 


o. , 

o., 

0.5, 

0 . , 

1. , 

o., 

1. , 

0.5, 

1. , 

1. , 

0.5, 

1-, 

0 . , 

1. , 

0. , 

0.5, 

0.5, 

0.5 

NODGl , 

INCR-X, -Y, -Z 




1, 

65, 

1, o, 





DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 
65, 


-1 

00, 

-0 

.995, 

-0 

99, 

-0 

98, 

-0 

97, 

-0 

96, 

-0 

95, 

-0 

935, 

-0 

92, 

-0 

905, 

-0 

89, 

-0 

87, 

-0 

85, 

-0 

83, 

-0 

81, 

-0 

785, 

-0 

76, 

-0 

73, 

-0 

70, 

-0 

66, 

-0 

62, 

-0 

575, 

-0 

53, 

-0 

48, 

-0 

43, 

-0 

38, 

-0 

33, 

-0 

28, 

-0 

23, 

-0 

175, 

-0 

12, 

-0 

06, 

0 

00, 

0 

06, 

0 

12, 

0 

175, 

0 

23, 

0 

28, 

0 

33, 

0 

38, 

0 

43, 

0 

48, 

0 

53, 

0 

575,. 

0 

62, 

0 

66 , 

0 

70, 

0 

73, 

0 

76, 

0 

785, 

0 

81, 

0 

83, 

0 

85, 

0 

87, 

0 

89, 

0 

905, 

0 

92, 

0 

935, 

0 

95, 

0 

96, 

0 

97, 

0 

98, 

0 

99, 

0 

995, 

1 

00, 















> 

-1 

00, 

-0 

995, 

-0 

99, 

-0 

98, 

-0 

97, 

-0. 

96, 

-0 

95, 

-0 

935, 

-0 

92, 

-0 

905, 

-0 

89, 

-0 

87, 

-0 

85, 

-0 

83, 

-0 

81, 

-0 

785, 

-0 

76, 

-0 

73, 

-0 

70, 

-0 

66, 

-0 

62, 

-0. 

575, 

-0 

53, 

-0 

48, 

-0 

43, 

-0 

38, 

-0 

33, 

-0 

28, 

-0 

23, 

-0. 

175, 

-0 

12, 

-0 

06, 

0 

oo, 

0 

06, 

0 

12, 

0 

175, 

0 

23 

0. 

28, 

0 

33, 

0 

38, 

0 

43, 

0 

48, 

0 

53, 

0 

575, 

0 

62, 

0 

66, 

0 

70, 

0 

73, 

0 

76, 

0 

785, 

0 

81, 

0 

83, 

0 

85, 

0 

87, 

0 

89, 

0 

905, 

0 

92, 

0 

935, 

0 

95, 

0 

96, 

0 

97, 

0. 

98, 

0 

99, 

0 

995 , 

1 

00, 















0 

i 
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ELEMENT CONNECTIVITY DATA FOR THE GLOBAL DOMAIN 

NUMBER OF BLOCKS (NBLOC) 

1 , 

ELEMENT NO. AND NODE CONNECTIVITY (IEL1 .NODES (IEL1) ) 

1. 1, 66, 131, 132, 133, 68, 3, 2, 67, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 

32, 32, 130,130,130, 130,130,130, 130,130,130, 

32, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 

IA01 --- 

INITIAL GUESS FOR U (NREC) 

0 , 

DBC FOR U 


4, 

1, 

4161, 

65, 

0. 

4161, 

4224, 

1, 

0. 

65, 

4225. 

65, 

1. 

1, 

64, 

1, 

0. 


IA02 

INITIAL GUESS FOR V (NREC) 
0 , 

DBC FOR V 
4, 


1, 

4161, 

65, 

0. , 

4161, 

4224, 

1, 

0., 

65, 

4225, 

65, 

o. , 

1, 

64, 

1, 

o., 

IA04 - - ( PBCDAT , 

IPNODl , 

IPNOD2) 

— 

0. , 

2017, 

2081, 



END OF INPUT DATA 

****PROCESSOR FOR NAVIER- STOKES EQUATIONS **** 

****END OF RUN 

******************************** BOTTOM OF DATA **************************** 



A. 2. 2 Backward-Facing Step Flow 


************************************** TOP OF DATA ************************* 
--- LAMINAR BACKWARD -FACING STEP FLOW (STP5K) --- 
****INITIALIZE DIMENSIONED VARIABLES **** 

****PREPARE INPUT DATA **** 

DESCRIPTIVE DATA 

I FLOW, NDIM, NGAUS , MFRONF , 

6, 2, 3, 95, 

CNTL PARAMETERS - --- 


NNODE , NELEM, IAXSY, IPLOT, 

2631, 628, 0, 1, 

ELEMENT CONNECTIVITY DATA FOR THE GLOBAL DOMAIN ----- 

NUMBER OF BLOCKS (NBLOC) 

3 , 

NODE CONNECTIVITY DATA FOR IBLOC=l (I ELI., NODES) 

1, 1, 16, 31, 32, 33, 18, 3, 2, 1.7, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 

3, 7, 30,30,30, 30,30,30, 30,30,30, 

7, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

1 , 0 , 0 , 0 , 0 , 0 ., 0 , 0 , 0 , 0 , 

. NODE CONNECTIVITY DATA FOR IBLOC=2 (IEL1, NODES) 

22, 91, 106, 137, 138, 139, 108, 93, 92, 107, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


0 , 


1, 

o, 


0, 0, 0, 0, 0, 0, 

0, 0, 0, 

7, 

1, 


2, 2, 2, 2, 2, 2, 

2, 2, 2, 

1, 

o. 


0, 0, 0, 0, 0, 0, 

0, 0, 0, 

NODE CONNECTIVITY DATA 

FOR 

IBLOG=3 (IEL1, NODES) 


29, 121, 

152, 183, 

184, 

185, 154, 123, 122, 

153, 

NEL(KDIM) , 

INCREL(KDIM) , 

INCNOD(KDIM) 


40, 

15, 


62,62,62, 62,62,62, 

62,62,62, 

15, 

1, 


2, 2, 2, 2, 2, 2, 

2, 2, 2, 

1, 

o, 


0, 0, 0, 0, 0, 0, 

0, 0, 0, 

NODE COORDINATE 

DATA 





NUMBER OF BLOCKS (NBLOC) 

2 , 

GRID GENERATION METHOD FOR IBLOC=l (METHOD) 
1, 

NODG1 , INCRX , INCRY , INCRZ, 


1 , 


15, 1, 


0, 


NDAT , GRID COORDINDATE DATA 


8, 

-0.0147, 

-0.0049, 

-0.01274, 

-0.00294, 

-0.01078, 

-0.00147, 

. -0.00882, 

-0.00686, 

15, 

0.0049, 

0.005145, 

0 .00539, 

0.0057085, 

0.006027, 


0.0064925, 

0.006958, 

0.0075, 

0.008042, 

0.0085075, 

1, 

0.008973, 
0 . , 

0.0092915, 

0.00961, 

0.009855, 

0.0101, 


GRID GENERATION METHOD FOR IBLOC=2 (METHOD) 

1 , 

NODG1 , INCRX, INCRY, INCRZ, 

121, 31, 1, 0, 

NDAT, GRID COORDINATE DATA 

81, 0., 0.00049, 0.00098, 0.00196, 

0.00441, 0.00588, 0.00784, 0.0098, 

0.014.7, 0.01715, 0.0196, 0.02205, 


0.00294, 

0.01225, 

0.0245, 



0.02695, 0.0294, 0.03185, 0.0343, 0.03675, 

0.0392, 0.04165, 0.0441, 0.04665, 0.049, 

0.05145, 0.0539, 0.05635, 0.0588, 0.06125, 

0.0637, 0.06615, 0.0686, 0.07105, 0.0735, 

0.07595, 0.0784, 0.08085, 0.0833, 0.08575, 

0.0882, 0.09065, 0.0931, 0.09555, 0.098, 

0.10045, 0.1029, 0.10535, 0.1078, 0.11025, 

0.1127, 0.11515, 0.1176, 0.12005, 0.1225, 

0.12495, 0.1274, 0.12985, 0.1323, 0.13475, 

0.1372, 0.13965, 0.1421, 0.14455, 0.147, 

0.14994, 0.15288, 0.15631, 0.15974, 0.16366, 

0.16758, 0.17199, 0.1764, 0.1813, 0.1862, 

0.19159, 0.19698, 0.202615, 0.20825, 0.214375, 

0.2205, 

31, 0., 0.000196, 0.000392, 0.0006615,0.000931, 

0.001274, 0.001617, 0.0020335,0.00245, 0.0028665, 

0.003283, 0.003626, 0.003969, 0.0042385,0.004508, 
0.004704, 0.0049, 0.005145, 0.00539, 0.0057085, 

0.006027, 0.0064925,0.006958, 0.0075, 0.008042, 

0.0085075,0.008973, 0.0092915,0.00961, 0.009855, 

0 . 0101 , 

1 , 0 . , 

MATERIAL PROPERTY OF FLUID --- (RE=500) - 

VISCY, DENSY, BFX(l-2), 

0.000016986, 1.225, 0. , 0. , 

ITERATION PARAMETERS - 

MAXIT, RELAX(l-lO), CNVCF(l-lO) , 

100 , 

0 . 8 , 0 . 8 , 1 ., 1 ., 1 ., 

1 ., 1 ., 1 ., 1 ., 1 ., 

l.E-4, l.E-4, l.E-4, l.E-4, l.E-4, 

l.E-4, l.E-4, l.E-4, l.E-4, l.E-4, 

IA01 

INITIAL GUESS FOR U-VELOCITY (NREC) 

0 , 

DBC FOR U 

20 , 


1, 

1, 

1, 

o.. 

2, 

2, 

1, 

0.1796, 

3, 

3, 

1, 

0.3414, 

4, 


1, 

0.5252, 

5, 

5, 

1, 

0.6790, 

6, 

6, 

1, 

0.8498, 

7, 

7, 

1, 

0.9565, 

8, 

8, 

1, 

1.0000, 

9, 

9, 

1, 

0.9565, 

10, 

10, 

1, 

0.8498, 

11, 

11, 

1, 

0.6790, 

12, 

12, 

1, 

0.5252, 

13, 

13, 

1, 

0.3414, 

14, 

14, 

1, 

0.1796, 

15, 

15, 

1, 

0. , 

16, 

106, 

15, 

o. , 

121, 

137, 

1, 

o.. 
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121,2601, 31, 0., 

30, 120, 15, 0. , 

151,2631, 31, 0., 

I AO 2 - 

INITIAL GUESS FOR V-VELOCITY (NREC) 


0, 



DBC FOR V 



6, 



1, 15, 

1, 

0 

16, 106, 

15, 

0 

121, 137, 

1, 

0 

121,2601, 

31, 

0 

30, 120, 

15, 

0 

151,2631, 

31, 

0 


IA04 ( PBCDAT , IPN0DE(l-2)) 

0., 0, 2617, 

END OF INPUT DATA 

****PROCESSOR FOR NAVIER- STOKES EQUATIONS **** 

****END OF RUN **** 

******************************** BOTTOM OF DATA ************************** 
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A. 2. 3 Flow Through a Nest of Cylinders 


********************************* TOP OF DATA ********************************* 

STEADY FLOW THROUGH A NEST OF CYLINDERS -(CNEST)- 

****INITIALIZE DIMENSIONED VARIABLES**** 

****PREPARE INPUT DATA **** 

DESCRIPTIVE DATA - --' - 

I FLOW , NDIM, NGAUS , MFRONF , 

6, 2, 3, 115, 

CNTL PARAMETERS 

NNODE , NELEM, IAXSY, I PLOT 
4369, 1024, 0,1, 

ELEMENT CONNECTIVITY DATA FOR THE GLOBAL DOMAIN 

NUMBER OF BLOCKS (NBLOC) 

21 , 

IBLOC=T , ( IEL1 , NODES(l-NPE) ) 

1, 1, 18, 35, 36, 37, 20, 3, 2, 19, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


20, 

8, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

8, 

1, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

1, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 


IBLOC=2 , ( IEL1 , NODES (1-NPE)) 

161, 681,698,715, 716,717,700, 683,682,699, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


8, 

8, 

34, 

34, 34, 

34, 34, 

34, 

34, 

34, 

34, 

8, 

1, 

2, 

2, 

2, 

2, 2, 

2, 

2, 

2, 

2, 

1, 

o, 

o, 

o, 

o, 

0, o, 

0, 

0, 

o, 

0, 

IBLOC=3 

-1, 

(IEL1 , 

NODES (1 

-NPE) ) 






225, 697,714,731, 972,989,988, 987,970,971, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


8, 

1, 

34, 

34, 

34, 

2, 2, 2, 

2, 

2, 

2, 

1, 

o, 

o, 

o, 

0, 

0, 0, 0, 

0, 

0, 

o, 

1, 

o, 

o, 

o, 

o, 

0, 0, 0, 

0, 

0, 

o, 

IBLOC=3 

-2, 

(IEL1 , 

NODES ( 1 ■ 

-NPE)) 




233, 

987,988,989, 

1006 , 1023 , 1022 , 

, 1021,1004,1005, 

NO. 

OF 

ELEMENTS (NEL, INCREL, INCNOD) 



8, 

1, 

2, 

2, 

2, 

2, 2, 2, 

2, 

2, 

2, 

7, 

8, 

34, 

34, 

34, 

34, 34, 34, 

34, 

34, 

34, 

1, 

o, 

o, 

0, 

0, 

0, 0, 0, 

0, 

0, 

o, 


IBLOC=4, (IEL1 , NODES (1-NPE)) 


289, 1225,1226,1227, 1244,1261,1260, 1259,1242,1243, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 

8 , 1 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 

8, 8, 34, 34, 34, 34, 34, 34, 34, 34, 34, 

1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 

IBLOC=5-l, (IEL1, NODES (1-NPE)) 

353, 1225,1242,1259, 1546,1547,1531, 1515,1514,1530, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


8, 

8, 

34, 

34, 

34, 

32, 

32, 

32, 

32, 

32, 

32, 

1, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

1, 

0, 

0, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 


IBLOC=5-2 , (IEL1, NODES(l-NPE) ) 

354, 1515,1531,1547, 1548,1549,1533, 1517,1516,1532, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 

8, 8, 32, 32, 32, 32, 32, 32, 32, 32, 32, 

7 , 1 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 
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1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 

IBLOC=6-l, IEL1 , NODES (1-NPE)) 

417, 1497,1786,1803, 1804,1805,1788, 1771,1770,1787, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1, 

o, 

o, 

0, 0, 

o, 

o, 

o, 

0, 

0, 

0, 

1. 

o, 

o, 

0, 0, 

o, 

o, 

o, 

0, 

0, 

0, 

1. 

o, 

o, 

0, 0, 

o, 

o, 

o, 

0, 

0, 

0, 

IBLOC=6 ■ 

■2, 

I ELI, 

NODES (1- 

NPE)) 







425, 1803,1820,1837, 1838,1839,1822, 1805,1804,1821, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


7, 

8, 

34, 

34, 

34, 

34, 

34, 34, 

34, 

34, 

34, 

1, 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

i; 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 


IBLOC=6 - 3 , I ELI , NODES (1-NPE) ) 

418, 1771,1788,1805, 1806 , 1807 , 1790 , 1773,1772,1789, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


8, 8, 

34, 

34, 34, 

34, 34, 

34, 

34, 

34, 

34, 

7, 1, 

2, 

2, 

2, 

2, 2, 

2, 

2, 

2, 

2, 

1, o, 

o, 

0, 

o, 

0, o, 

0, 

0, 

0, 

0, 

IBLOC=7 - 1 , 

(IEL1 , 

NODES (1- 

-NPE)) 






481, 1497,2058,2059, 2075 , 2091 , 2090 , 1803 , 1786 , 2074 , 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

1 . 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 


IBLOC=7-2 , (IEL1 , NODES (1-NPE) ) 

489, 1803,2090,2091, 2107,2123,2122, 1837,1820,2106, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

0, 

8, 

34, 

32, 

32, 

32, 

32, 

32, 

34, 

34, 

32, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

0, 

0, 

0, 


IBLOC=7 - 3 , (I ELI, NODES (1 -NPE) ) 

482, 2059,2060,2061, 2077,2093,2092, 2091,2075,2076, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


7, 

1, 

2, 

2, 

2, 

2, 2, 

2, 

2, 

2, 

2, 

8, 

8, 

32, 

32, 32, 

32, 32, 

32, 

32, 

32, 

32, 

1, 

o, 

o, 

0, 

o, 

0, o, 

0, 

0, 

0, 

0, 

IBLOC=8- 

■1, 

(IEL1 , 

NODES (1 

-NPE)) 






545, 2041,2314,2315, 2332,2349,2348, 2347,2330,2331, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

1, 

o. 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

1, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 


IBLOC=8-2 , (IEL1 , NODES(l-NPE) ) 

553, 2347,2348,2349, 2366,2383,2382, 2381,2364,2365, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1, 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 0 , 0 , 

0 , 

7, 

8, 

34, 

34, 

34, 

34, 

34, 

34, 34, 34, 

34, 

1, 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 0 , 0 , 

0 , 


IBLOC=8 - 3 , (IEL1 , NODES (1-NPE) ) 

546, 2315,2316,2317, 2334,2351,2350, 2349,2332,2333, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 

7, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

8, 8, 34, 34, 34, 34, 34, 34, 34, 34, 34, 
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1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 

IBLOC=9 - 1 , (IEL1 , NODES (1-NPE)) 

609, 2601,2602,2619, 2620,2621,2604, 2599,2600,2603, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1 , 

1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

8 , 

1 , 

- 2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

- 2 , 

- 2 , 

2 , 

1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 


IBLOC=9-2 , (IEL1 , NODES (1-NPE) ) 

617, 2619,2636,2653, 2654,2655,2638, 2621,2620,2637, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


7, 

8, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

8, 

1, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

1, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 


IBLOC=10-1, (IEL1 , NODES (1-NPE)) 

673, 2873,2874,2891, 2892,2893,2876, 2839,2856,2875, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1, o, 

0, 

0, 0, 

o, o, 

o, 

0, 

0, 

0, 

7, 1, 

-34, 

2, 2, 

2, 2, 

2, 

-34, 

-34, 

2, 

1. o, 

0, 

0, 0, 

o, o, 

O', 

0, 

0, 

0, 

IBLOC=10-2 , 

(IEL1 , 

NODES (1-NPE)) 






680, 2635,2888,2905, 2906,2907,2890, 2585,2618,2889, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


1, 

o, 

0, 

o, 

o, 

o, 

0, 

0, 

0, 

0, 

0, 

1, 

o, 

0, 

0, 

o, 

o, 

o, 

0, 

0, 

0, 

0, 

1, 

o, 

0, 

0, 

o, 

o, 

0, 

0, 

0, 

0, 

0, 

IBLOC=10 

-3, 

(IEL1 , 

NODES (1 

-NPE) ) 






681, 2891,2908,2925, 2926,2927,2910, 2893,2892,2909, 

NO. OF ELEMENTS (NEL, INCREL, INCNOD) 


43, 

8, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

34, 

8, 

1, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

2, 

1, 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 


NODE COORDINATE DATA 

NUMBER OF BLOCKS (NBLOC) 

10 , 

GRID GENERATION METHOD FOR IBLOC=l 

2 , 

NODE COORDINATE DATA ( (XNOD(KPE,KDIM) ,KDIM=1 ,NDIM) ,KPE=1 ,NPE) 


o., 

o., 

1.25, 

0. , 

2.5, 

o., 

2.5, 

0.5,. 

2.5, 

1., 

1.25, 

1. , 

o., 

1., 

0., 

0.5, 

1.25, 

0.5, 




NODG1 , INCR-X, -Y, -Z 
1, 17, 1, 0, 

DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 
41. 


-1 

» 

-0 

95, 

-0 

9, 

-0 

85, 

-0 

8, 

-0 

75, 

-0 

7, 

-0 

65, 

-0 

6, 

-0 

55, 

-0 

50, 

-0 

45, 

-0 

40, 

-0 

35, 

-0 

30, 

-0 

25, 

-0 

20, 

-0 

15, 

-0 

10, 

-0 

05, 

0 

l 

0 

05, 

0 

10, 

0 

15, 

0 

20, 

0 

25, 

0 

30, 

0 

35, 

0 

40, 

0 

45, 

0 

50, 

0 

55, 

0 

60, 

0 

65, 

0 

70, 

0 

75, 

0 

80, 

0 

85, 

0 

90, 

0 

95, 
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1 . 0 , 


17 , 

-1., -0.875, -0.75, 

-0.375, -0.250, -0.125, 

0.25, 0.375, 0.50, 

0.875, 1.0, 

1 , 

0., 

GRID GENERATION METHOD FOR IBLOC=2 

2 , 

NODE COORDINATE DATA ((XNOD(KPE 
2.5, 0., 

3.0, 0., 

3.14645, 0.35355, 

2.5, 1., 

2.76903, 0.34567 

NODG1 , INCR-X, -Y, -Z 
681, 17, 1, 0, 


-0.625, -0.50, 

0.0, 0.125, 

0.625, 0.750, 


KDIM) , KDIM=1 , NDIM) , KPE=1 , NPE) 
2.75, 0., 

3.03806, 0.19134, 

2.82323, 0.67678, 

2.5, 0.5, 


DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT, DELX- ARRAY) 
17, 


-1., 

-0.8, 

-0.6, 

-0.42, 

-0.24, 

-0.08, 

0.08, 

0.22, 

0.36, 

0.48, 

0.60, 

0.96, 

17, 

0.70 

1., 

0.80, 

0.86, 

0.92, 

-1., 

-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

o.o, 

0.125,. 

0.25, 

0.875, 

1, 

o., 

0.375, 

1.0, 

0.50, 

0.625, 

0.750, 


u - i 

GRID GENERATION METHOD FOR IBLOC=3 

2 , 


NODE COORDINATE 

DATA ( (XNOD (KPE , KDIM) , KDIM=1 , 

NDIM) , KPE=1 ,NPE) 

2-5, 

1. , 

2.82323, 

0.67678, 

3.14645, 

0.35355, 

3.30866, 

0.46194, 

3.5, 

0.5, 

3.5, 

0.75, 

3.5, 

1., 

3. , 

1., 

3.15433, 

0.73097, 



NODG1 , INCR-X, 

-Y,-Z 



970, 1, 

17, 0, 




DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 


17, 


-1., 

-0.8, 

-0.6, 

-0.42, 

-0.24, 

-0.08, 

0.08, 

0.22, 

0.36, 

0.48, 

0.60, 

0.70 

0.80, 

0.86, 

0.92, 

0.96, 

1., 




16, 

-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

0 . 0 , 

0.125, 

0.25, 

0.375, 

0.50, 

0.625, 

0.750, 


0.875, 1.0, 

1 , 

0 ., 
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GRID GENERATION METHOD FOR IBL0C=4 


2 , 

NODE COORDINATE DATA ( (XNOD(KPE , KDIM) ,KDIM=1 , NDIM) ,KPE=1 , NPE) 


3.5, 

3.5, 

4., 

4.5, 

3.926775, 


1 . . 

0.5, 

0 ., 

0 ., 

0.426775, 


3.5, 

3.85355, 

4.25, 

4., 


0.75, 
0.35355, 
0 . , 

0.5, 


NODG1 , INCR-X, -Y, -Z 
1225, 1, 17, 0, 

DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 


17, 


-1.. 

-0.8, 

-0.6, 

-0.42, 

-0.24, 

-0.08, 

0.08, 

0.22, 

0.36, 

0.48, 

0.60, 

0.96, 

17, 

0. 70 

1. , 

0.80, 

0.86, 

0.92, 

-1.0, 

-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

0.0, 

0.125, 

0.25, 

0.875, 

1. 

0 . , 

0.375, 

1.0, 

0.50, 

0.625, 

0.750, 


v • i 

GRID GENERATION METHOD FOR IBLOC=5 
2 , 


NODE COORDINATE DATA ( (XNOD(KPE , KDIM) ,KDIM=1 , NDIM) , KPE=1 , NPE) 


3.5, 

1., 


4., 

0.5, 

4.5, 

0., 


4.5, 

0.25, 

4.5, 

0.5, 


4.14645, 

0.64645, 

4. , 

1., 


3.75, 

1. , 

4.07322, 

0.573225, 



NODG1 , INCR 

-X.-Y.-Z 




1514, 

16, 1, 0, 




DISCRETIZATION OF THE 

COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 

17, 





-i.o, 

-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

0.0, 

0.125, 

0.25, 

0.375, 

0.50, 

0.625, 

0.750, 

0.875, 

1.0, 




16, 

-0.8, 

-0.6, 

-0.42, 

-0.24, 

-0.08, 

0.08, 

0.22, 

0.36, 

0.48, 

0.60, 

0.70 

0.80, 

0.86, 

0.92, 

0.96, 

1. . 




1, 





0., 





GRID GENERATION 

METHOD FOR IBLOC=6 



. 

NODE COORDINATE DATA 

( (XNOD(KPE 

, KDIM) ,KDIM=1, 

NDIM) ,KPE=1, NPE) 

4.5, 

o., 


5. , 

0.5, 

5.5, 

1., 


5.25, 

1., 

5. , 

1-, 


4.85355, 

0.64645, 

4.5, 

0.5, 


4.5, 

0.25, 

4.926775, 

0.573225, 





NODGl, INGR-X, -Y, -Z 
1786, 17, 1, 0, 

DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT , DELX -ARRAY) 
16, 



-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250,. 

-0.125, 

0.0, 

0.125, 

0.25, 

0 . 875 , 

0 . 375 , 
1.0, 

0.50, 

0 . 625 , 

0.750, 

17, 

-1.0, 

-0.8V 

-0.6, 

-0.42, 

-0.24, 

-0.08, 

0.08, 

0.22, 

0.36, 

0.48, 

0.60, 

0.96, 

0. 70 

1. . 

0.80, 

0.86, 

0.92, 

1, 

o. , 

GRID GENERATION 

METHOD 

FOR IBLOC=7 



2, 

NODE COORDINATE DATA 

( (XNOD (KPE , KDIM) , KDIM= 

1,NDIM) ,KPE= 

4.5, 

0 . , 


4.75, 

0 ., 

5. , 

o., 


5.14645, 

0.35355 

5.5, 

0.5, 


5.5, 

0.75, 

5.5, 

1 . , 


5. , 

0.5, 


5.073225, 0.426775, 

NODGl , INCR-X,,-Y, -Z 
2058, 1, 16, 0, 

DISCRETIZATION OF THE COMPUTATIONAL GRID (NDAT, DELX -ARRAY) 
16, 



-0.8, 

-0.6, 

-0.42, 

-0.24, 

-0 . 08 , 

0.08, 

0.22, 

0.36, 

0.48, 

0.60, 

0.96, 

17, 

O'. 70 
1. , 

0,80, 

0.86, 

0.92, 

-1., 

-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

0.0, 

0.125 

0.25, 

0.875, 

0.375, 

1.0, 

0.50, 

0.625, 

0.750 


1 , 

0 . 


GRID GENERATION METHOD FOR IBLOC=8 

O 



. 

NODE COORDINATE DATA 

((XNOD (KPE, KDIM) ,KDIM=1,NDIM) ,KPE=1,NPE) 

5.5, 

1. , 


5.5, 

0.75, 

5.5, 

0.5, 


5.69134, 

0.46194, 

5. 85355 , 

0.35355, 

6.176775, 

0.676775, 

6.5, 

1., 


6. , 

1. , 

5 . 84567 , 

0.73097, 



NODGl, INCR-X, 

,-Y,-Z 




2330, 1, 

17, 0, 




DISCRETIZATION 

1,7 

OF THE 

COMPUTATIONAL GRID (NDAT, DELX- ARRAY) 

1/ , 

-1., -0. 

,8, 

-0.6, 

-0.42, 

-0.24, 

-0.08, 0. 

.08, 

0.22, 

0.36, 

0.48, 

0.60, 0. 

70 

0.80, 

0.86, 

0 . 92 , 

0.96, 1. 

• > 
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16 , 


-0.875, -0.75, -0.625, -0.50, 

-0.375, -0.250, -0.125, 0.0, 0.125, 

0.25, 0.375, 0.50, 0.625, 0.750, 

0.875, 1.0, 

1 , 

0 . , 

GRID GENERATION METHOD FOR IBLOC=9 

2 , 


NODE COORDINATE DATA 

((XNOD(KPE.KDIM) .KDIM- 

l.NDIM) ,KPE=1,NPE) 

5.85355, 

0.35355, 

5.96194, 

0.19134, 

6., 

o., 


6.25, 

0 . , 

6.5, 

o., 


6.5, 

0.5, 

6.5, 

1., 


6.176775, 0.676775, 

6.23097, 

0.34567, 



NODG1 , INCR 

-X.-Y.-Z 




2602, 

17, 1, 0, 




DISCRETIZATION OF THE 

COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 

16, 






-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

0.0, 

0.125, 

0.25, 

0.375, 

0.50, 

0.625, 

0.750, 

0.875, 

1.0, 




17, 





-1.. 

-0.96, 

-0.92, 

-0.86, 

-0.80, 

-0.70, 

-0.60, 

-0.48, 

-0.36, 

-0.22, 

-0.08, 

0.08, 

0.24, 

0.42, 

0.6, 

0.8, 

1.0, 




1, 





0 ., 





GRID GENERATION 
o 

METHOD FOR IBL0C=10 



^ * 

NODE COORDINATE DATA 

( (XNOD (KPE , KDIM) , KDIM= 

d.NDIM) , KPE=1,NPE) 

6.5, 

o., 


23.75, 

0 . , 

41.0, 

0 ., 


41. , 

0.5, 

41.0, 

1., 


23.75, 

1., 

6.5, 

1., 


6.5, 

0.5, 

23.75, 

0.5, 




NODG1 , INCR 

-X.-Y.-Z 




2874, 

17, 1, 0, 




DISCRETIZATION OF THE 

COMPUTATIONAL GRID (NDAT , DELX- ARRAY) 

88, 





-0.997, 

-0.994, 

-0.991, 

-0.988, 

-0.985, 

-0.982, 

-0.979, 

-0.976, 

-0.973, 

-0.97, 

-0.967, 

-0.964, 

-0.961, 

-0.958, 

-0.954, 

-0.95, 

-0.945, 

-0.94, 

-0.934, 

-0.928, 

-0.921, 

-0.914, 

-0.906, 

-0.898, 

-0.889, 

-0.88, 

-0.87, 

-0.86, 

-0.848, 

-0.836, 

-0.823, 

-0.81, 

-0.794, 

-0.778, 

-0.76, 

-0.742, 

-0.722, 

-0.702, 

-0.68, 

-0.658, 

-0.635, 

-0.612, 

-0.588, 

-0.564, 

-0.54, 

-0.516, 

-0.492, 

-0.468, 

-0.443, 

-0.418, 

-0.393, 

-0.368, 

-0.34, 

-0.312, 

-0.282, 

-0.252, 

-0.221, 

-0.19, 

-0.155, 

-0.12, 
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-0.08, 

-0.04, 

0.0, 

0.04, 

0.08, 

0.12, 

0.16, 

0.20, 

0.24, 

0.28, . 

0.32, 

0.36, 

0.40, 

0.44, 

0.48, 

0.52, 

0.56, 

0.60, 

0.64, 

0.68, 

0.72, 

0.92, 

0.76, 

0.96, 

0.80, 

1., 

0.84, 

0.88, 

-1.0, 

-0.875, 

-0.75, 

-0.625, 

-0.50, 

-0.375, 

-0.250, 

-0.125, 

0.0, 

0.125, 

0.25, 

0.875, 

o.. 

0.375, 

1.0, 

0.50, 

0.625, 

0.750, 


MATERIAL PROPERTY OF THE FLUID 

VISCY (RE=40) , DENSY, BFX(1,2) 
0.030625, 1.225, 0., 0., 


ITERATION PARAMETERS 


MAXIT, RELAX(l-lO) , CNVCF(l-lO) 

50, 

0 . 8 , 0 . 8 , 1 ., 1 ., 1 ., 

1 ., 1 ., 1 ., 1 ., 1 ., 

l.E-4, l.E-4, l.E-4, l.E-4, l.E-4, 

l.E-4, l.E-4, l.E-4, l.E-4, l.E-4, 

1A01 - 


NREC FOR INITIAL GUESS 
0 , 

DBC FOR U 


9, 


1, 

17, 

1, 

1., 

953, 

969, 

1, 

0., 

986, 

1241, 

17, 

0., 

1258, 

1513, 

17, 

o., 

1529, 

1785, 

16, 

o., 

1802, 

2057, 

17, 

o., 

2073, 

2329, 

16, 

o. , 

2346, 

2601, 

17, 

o. , 

2602, 

2857, 

17, 

o. , 


I AO 2 - - 

NREC FOR INITIAL GUESS 

0 , 

DBC FOR V 


21 , 


1, 

17, 

1, 

o.. 

1, 

681, 

17, 

o.. 

17, 

697, 

17, 

o. , 

698, 

953, 

17, 

o. , 

953, 

969, 

1, 

o. , 

970, 

1225, 

17, 

0., 

986, 

1241, 

17, 

o. , 

1258, 

1513, 

17, 

o. , 

1497, 

1513, 

1, 

0. , 

1514, 

1529, 

1, 

o. , 

1529, 

1785, 

16, 

0., 

1802, 

2057, 

17, 

o. , 



2041, 

2057, 

1. 

o. , 

2058, 

2073, 

1, 

o., 

2073, 

2329, 

16, 

0. , 

2330, 

2585, 

17, 

o. , 

2346, 

2601, 

17, 

o. , 

2602, 

2857, 

17, 

o., 

2857, 

2873, 

1. 

o., 

2874, 

4353, 

17, 

o., 

2890, 

4369, 

17, 

o., 


IA04 - - ( PBCDAT , IPN0D1, IPNOD2) 

0., -953, 953, 

END OF INPUT DATA 

**** PROCESSOR FOR NAVIER- STOKES EQUATIONS **** 

****end **** 

******************************** BOTTOM OF DATA ******************************** 
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A. 2.4 Channel Flow with an Internal Blockage 


************************************* TOP OF DATA ********************************* 
--- CHANNEL FLOW WITH BLOCKAGE --- BL0CK5 --- 

****initialize dimensioned variables **** 

****PREPARE input data **** 

DESCRIPTIVE DATA - - - - 

I FLOW, NDIM, NGAUS , HFRONF , 

6, 2, 3, 95, 

CNTL PARAMETERS --- ----- 

NNODE , NELEM, IAXSY, I PLOT, 

2085 490, O', 1, 

ELEMENT CONNECTIVITY DATA FOR THE GLOBAL DOMAIN 

NUMBER OF BLOCKS (NBLOC) 

13 


NODE CONNECTIVITY DATA FOR IBLOC=l AND 2 (I ELI, NODES) - B1 
1, 1, 26, 51, 52, 53, 28, 3, 2, 27, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


10, 

12, 

50,50,50, 

50,50,50, 

50,50,50, 

12, 

1, 

2, 2, 2, 

2, 2, 2, 

2, 2, 2, 

1, 

0, 

0, 0, 0, 

0, 0, 0, 

0, 0, 0, 

NODE CONNECTIVITY DATA FOR 

IBLOC=3-l 

(IEL1, NODES) - B2 

121, 501, 
NEL(KDIM) 

526, 541, 542, 
, INCREL(KDIM) , 

543, 528, 503, 502, 
INCNOD(KDIM) 

527, 

1. 

o, 

o, 0, 0, 

0, 0, 0, 

0, 0, 0, 

7, 

1, 

2, 2, 2, 

2, 2, 2, 

2, 2, 2, 

1, 

o, 

0, 0, 0, 

0, 0, 0, 

0, 0, 0, 


NODE CONNECTIVITY DATA FOR IBLOC=3-2 (IEL1, NODES) - B3 
128, 541, 556, 571, 572, 573, 558, 543, 542, 557, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

7, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 

NODE CONNECTIVITY DATA FOR IBLOC=4-l (IEL1, NODES) - B4 
135, 585,586,601, 602,603,588, 555,570,587, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


o, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 


NODE CONNECTIVITY DATA FOR IBLOC=4-2 (IEL1, NODES) - B5 
136, 555,588,603, 604,605,590, 515,540,589, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

1 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 


NODE CONNECTIVITY DATA FOR IBLOC=5-l (I ELI, NODES) - B6 
137, 515,590,605, 606,607,592, 517,516,591, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


1, 

0, 

0, 

0, 0, 

0, 

0, 

0, 

0, 0, 

o, 

5, 

1, 

2, 

2, 2, 

2, 

2, 

2, 

2, 2, 

2, 

1. 

0, 

0, 

0, o, 

0, 

0, 

0, 

0, 0, 

o, 


NODE CONNECTIVITY DATA FOR IBLOC=4-3 AND 5-2 (IEL1, NODES) - B7 
142, 601,616,631, 632,633,618, 603,602,617, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 

5, 7, 30,30,30, 30,30,30, 30,30,30, 

7 , 1 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 , 
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1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 
NODE CONNECTIVITY DATA FOR IBLOC=6-l (IEL1, NODES) - B8 
177, 766,780,794, 795,796,782, 768,767,781, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


2, 

7, 

28,28,28, 

28,28,28, 

28,28,28, 

6, 

1, 

2, 2, 2, 

2, 2, 2, 

2, 2, 2, 

1, 

0, 

0, 0, 0, 

0, 0, 0, 

0, 0, 0, 


NODE CONNECTIVITY DATA FOR IBLOC=6-2 (IEL1, NODES) - B9 
183, 778,792,806, 807,753,752, 751,779,793, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


2, 

7, 

28,28,28, 

28, 

2, 

2, 

2, 

28, 

28, 

1, 

0 , 

0 , 0 , 0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

1, 

0 , 

0 , 0 , 0 , 

0 , 

0 , 

0 , 

0 , 

0 , 

0 , 


NODE CONNECTIVITY DATA FOR IBLOC=7-l (IEL1, NODES) - BIO 
191, 822,836,861, 862,863,838, 824,823,837, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


1, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

6, 

1, 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

1, 

0, 

0, 

0, 

0, 

o, 

0, 

0, 

0, 

o, 

0, 


NODE CONNECTIVITY DATA FOR IBLOC=7-2 (IEL1, NODES) - Bll 
197, 834,848,873, 874,875,850, 755,835,849, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

0, 


NODE CONNECTIVITY DATA FOR IBLOC=8-l (IEL1, NODES) - B12 
198, 755,850,875, 876,877,852, 757,756,851, 

NEL(KDIM), INCREL(KDIM) , INCNOD(KDIM) 


0, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

0, 

0, 

1, 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

2 , 

0, 

0, 

0, 

0, 

0, 

0, 

0, 

o, 

0, 

o. 


NODE CONNECTIVITY DATA FOR IBLOC=7-3 AND 8-2 (IEL1, NODES) - B13 
203, 861,886,911, 912,913,888, 863,862,887, 

NEL(KDIM) , INCREL(KDIM) , INCNOD(KDIM) 


24, 

12, 

50,50,50, 

50,50,50, 

50,50,50, 

12. 

1, 

2, 2, 2, 

2, 2, 2, 

2, 2, 2, 

1, 

0, 

0, 0, 0, 

0, 0, 0, 

0, 0, 0, 


NODE COORDINATE DATA 


NUMBER OF BLOCKS (NBLOC) 

8 , 

1. GRID GENERATION METHOD FOR IBLOC=l (METHOD) 

2 , 


GLOBAL NODE COORD. DATA (XNOD(KPE , KDIM) , KDIM=1 , NDIM) , KPE=1 , NPE) 


0. , 0. , 3.375, 0. , 

6.75,0.625, 6.75, 1.25, 

0., 1.25, 0., 0.625, 


6.75, 0., 

3.375, 1.25, 
3.375, 0.625, 


NODG1 , INCRX , INCRY , INCRZ , 


1, 

25, 

1, 


0, 





NDAT , 

GRID 

COORDINDATE 

DATA 




21, 

-1. 

» 

-0. 

.852, 

-0, 

.704, 

-0. 

.555, 


-0. 

258, 

-0. 

.110, 

0. 

.038, 

0. 

,186, 


0. 

438, 

0. 

,534, 

0. 

.630, 

0. 

,704, 


0. 

830, 

0. 

,882, 

0. 

,919, 

0. 

,956, 


1. 

0, 








-0.406, 

0.312, 

0.778, 

0.978, 
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15, -1., -0.92, -0.84, -0.56, -0.28, 

-0.04, 0.20, 0.32, 0.44, 0.52, 

0.60, 0.68, 0.76, 0.88, 1.0, 

1 , 0 ., 

2. GRID GENERATION METHOD FOR IBLOC=2 (METHOD) 

2 , 

GLOBAL NODE COORD. DATA (XNOD(KPE , KDIM) ,KDIM=1 , NDIM) ,KPE=1,NPE) 
0., 1.25, 3.375,1.25, 6.75, 1.25, 

6.75,1.875, 6.75, 2.5, 3.375,2.5, 

0., 2.5, 0., 1.875, 3.375,1.875, 

NODG1 , INCRX , INCRY , INCRZ , 


15, 

25, 1, 

0 , 




NDAT, 

GRID COORDINDATE 

DATA 



21, 

-1., 

-0.852, 

-0.704, 

-0.555, 

-0.406, 


-0.258, 

-0.110, 

0.038, 

0.186, 

0.312, 


0.438, 

0.534, 

0.630, 

0.704, 

0.778, 


0.830, 

1.0, 

0.882, 

0.919, 

0.956, 

0.978, 

11, 

-1-, 

-0.88, 

-0.76, 

-0.52, 

-0.28, 


0.0, 

0.28, 

0.56, 

0.84, 

0.92, 


1 . 0 , 

1 , 0 ., 

3. GRID GENERATION METHOD FOR IBLOC=3 (METHOD) 

2 , 

GLOBAL NODE COORD. DATA (XNOD (KPE, KDIM) ,KDIM=1 , NDIM) ,KPE=1 ,NPE) 

6.75.0. , 6.875,0., 7., 0., 

7., 0.5, 7., 1., 6.875,1.125, 

6.75,1.25, 6.75, 0.625, 6.875,0.5625, 

NODG1 , INCRX , INCRY , INCRZ , 


526, 

NDAT, 

15, 1, 0, 

GRID COORDINDATE 

DATA 



4, 


-0.4, 

0.2, 

0.6, 

i.o, 

15, 

-1. , 

-0.92, 

-0.84, 

-0.56, 

-0.28, 


-0.04, 

0.20, 

0.32, 

0.44, 

0.52, 

1, 

0.60, 

o., 

0.68, 

0.76, 

0.88, 

1.0, 

GRID GENERATION 

METHOD 

FOR IBL0C=4 

(METHOD) 



2 , 


GLOBAL NODE COORD. DATA (XNOD (KPE, KDIM) ,KDIM=1 , NDIM) ,KPE=1 ,NPE) 


7. , 

1. , 

7.5, 1., 

8. , 

1., 

8.125, 

1.125, 

8.25, 1.25, 

7.5, 

1.25, 

6.75, 

1.25, 

6.875,1.125, 

7.5, 

1.125 


NODG1 , INCRX , INCRY , INCRZ , 

586, 15, 1, 0, 

NDAT, GRID COORDINDATE DATA 

12, -0.83, -0.67, -0.5, -0.33, 

-0.17, 0., 0.17, 0.33, 0.5, 

0.67, 0.83, 1., 

5, -1., -0.6, -0.2, 0.4, 1., 

1 , 0 ., 

5. GRID GENERATION METHOD FOR IBLOC=5 (METHOD) 

2 , 

GLOBAL NODE COORD. DATA (XNOD (KPE , KDIM) ,KDIM=1 , NDIM) ,KPE=1 ,NPE) 
6.75, 1.25, 7.5, 1.25, 8.25, 1.25, 
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8.25, 

1.875, 

8.25, 

2.5, 

7.5, 2.5, 

6.75, 

2.5, 

6.75, 

1.875, 

7.5, 1.875, 

NODG1 , INCRX , INCRY , INCRZ , 



590, 

15, 1, 

o, 



NDAT, GRID COORDINDATE DATA 


12, 


-0.83, 

-0.67, 

-0.5, -0.33, 


-0.17, 

o., 

0.17, 

0.33, 0.5, 


0.67, 

0.83, 

■ 1. , 


11, 

-1. , 

-0.88, 

-0.76, 

-0.52, -0.28, 


0.0, 

0.28, 

0.56, 

0.84, 0.92, 


1.0, 




1, 

o., 




GRID GENERATION METHOD FOR IBLOC=6 

(METHOD) 

2, 





GLOBAL 

NODE COORD. DATA 

(XNOD (KPE, KDIM) ,KDIM=1, NDIM) 

8., 

0. , 

8.125, 0 

8. 

25, 0., 

8.25, 

0.625, 

8.25, 1 

.25, 8. 

125, 1.125, 


8. , 1. , 8. , 0.5, 8.125, 0.5625, 

NODG1 , INCRX , INCRY , INCRZ , 

766, 14, 1, 0, 


NDAT , GRID COORDINDATE DATA 


5, 

-1. , 

-0.6, 

■0.2, 

0.4, 

i.o, 

14, 

-1. , 

-0.92, 

-0.84, 

-0.56, 

-0.28, 


-0.04, 

0.20, 

0.32, 

0.44, 

0.52, 


0.60, 

0.68, 

0.76, 

0.88, 


1, 

o., 





GRID GENERATION 

METHOD FOR 

IBLOC=7 

(METHOD) 



GLOBAL NODE COORD. DATA (XN0D(KPE , KDIM) , KDIM-1 ,NDIM) , KPE=1 , NPE) 
8.25, 0. , 19.125, 0. , 30. , 0. , 

30.0, 0.625, 30.0, 1.25, 19.125, 1.25, 

8.25, 1.25, 8.25, 0.625, 19.125, 0.625, 

NODG1 , INCRX , INCRY , INCRZ , 

836, 25, 1, 0, 

NDAT, GRID COORDINDATE DATA 


-0 

993, 

-0 

986, 

-0 

975, 

-0 

964, 

-0 

948, 

-0 

932, 

-0 

909, 

-0 

886, 

-0 

856, 

-0 

826, 

-0 

787, 

-0 

748, 

-0 

702, 

-0 

656, 

-0 

610, 

-0 

564, 

-0 

518, 

-0 

472, 

-0 

426, 

-0 

380, 

-0 

334, 

-0 

288, 

-0 

242, 

-0 

196, 

-0 

150, 

-0 

104, 

-0 

058, 

-0 

012, 

0 

034, 

0 

080, 

0 

126, 

0 

172, 

0 

218. 

0 

264, 

0 

310, 

0 

356, 

0 

402, 

0 

448, 

0 

494, 

0 

540, 

0 

586, 

0 

632, 

0 

678, 

0 

724, 

0 

770, 

0 

816, 

0 

862, 

0 

908, 

0 

954, 

1 

o, 

-1 

» 

-0 

92, 

-0 

84, 

-0 

56, 

-0 

28, 

-0 

04, 

0 

20, 

0 

32, 

0 

44, 

0 

52, 

0 

60, 

0 

68, 

0 

76, 

0 

88, 

1 

0 , 


1 , 0 ., 

8. GRID GENERATION METHOD FOR IBL0C=8 (METHOD) 

2 , 

GLOBAL NODE COORD. DATA (XNOD(KPE , KDIM) , KDIM=1 , NDIM) ,KPE=1 , NPE) 
8.25, 1.25, 19.125, 1.25, 30.0, 1.25, 

30.0, 1.875, 30.0, 2.5, 19.125, 2.5, 
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8.25, 2.5, 


8.25, 1.875, 19.125, 1.875, 


N0DG1 , INCRX , INCRY , INCRZ , 
850, 25, 1, 0, 

NDAT, GRID COORDINDATE DATA 


-0 

993, 

-0 

986, 

-0 

975, 

-0 

964, 

-0 

948, 

-0 

932, 

-0 

909, 

-0 

886, 

-0 

856, 

-0 

826, 

-0 

787, 

-0 

748, 

-0 

702, 

-0 

656, 

, -0 

610, 

-0 

564, 

-0 

518, 

-0 

472, 

-0 

426, 

-0 

380, 

-0 

334, 

-0 

288, 

-0 

242, 

-0 

196, 

-0 

150, 

-0 

104, 

-0 

058, 

-0 

012, 

0 

034, 

0 

080, 

0 

126, 

0 

172, 

0 

218, 

0 

264, 

0 

310, 

0 

356, 

0 

402, 

0 

448, 

0 

494, 

0 

540, 

0 

586, 

0 

632, 

0 

678, 

0 

724, 

0 

770, 

0 

816, 

0 

862, 

0 

908, 

0 

954, 

1 

0, 

-1 

> 

-0 

88, 

-0 

76, 

-0 

52, 

-0 

28, 

0 

0, 

0 

28, 

0 

56, 

0 

84, 

0 

92, 


1 . 0 , 

1 , 0 . , 

MATERIAL PROPERTY OF FLUID --- (RE=500) 


VISCY, DENSY, BFX(l-2) , 
0.006125, 1.225, 0. , 0. , 

ITERATION PARAMETERS --- 


IA01 


MAXIT, RELAX(l-lO), CNVCF(l-lO), 


50, 

0 . 8 , 0 . 8 , 

1-, I-, 

l.E-6, l.E-6, 

l.E-4, l.E-4, 


1., I-, I-, 

1., 1., I-, 

l.E-4, l.E-6, l.E-4, 

l.E-4, l.E-4, l.E-4, 


INITIAL GUESS FOR U-VELOCITY (NREC) 

0 , 

DBC FOR U 
35, 


1, 

1, 

1, 

o. , 

2, 

2, 

1, 

0.0784, 

3, 

3, 

1, 

0.1536, 

4, 

4, 

1, 

0.3916 

5, 

5, 

1, 

0.5904, 

6, 

6, 

1, 

0.7296, 

7, 

7, 

1, 

0.84, 

8, 

8, 

1, 

0.8844, 

9, 

9, 

1, 

0.9216, 

10, 

10, 

1, 

0.9424, 

11, 

11, 

1, 

0.96, 

12, 

12, 

1, 

0.9744, 

13, 

13, 

1, 

0.9856, 

14, 

14, 

1, 

0.9964, 

15, 

15, 

1, 

i.o, 

16, 

16, 

1, 

0.9964, 

17, 

17, 

1, 

0.9856, 

18, 

18, 

1, 

0.9424, 

19, 

19, 

1, 

0.8704, 

20, 

20, 

1, 

0.75, 

21, 

21, 

1, 

0.5904, 



22 , 

22 , 

1 , 

0 . 3916 , 

23 , 

23 , 

1 . 

0 . 1536 , 

24 , 

24 , 

1 , 

0 . 0 . 0784 , 

25 , 

25 , 

1 . 

o. , 

1 , 501 , 

25 , 

o. , 

25 , 525 , 

25 , 

o. , 

526 , 571 , 

15 , 

o. , 

571 , 585 , 

1 . 

o. , 

586 , 751 , 

15 , 

o. , 

600 , 765 , 

15 , 

o. , 

766 , 779 , 

1 . 

o. , 

780 , 822 , 

14 , 

o., 

836 , 2061 , 

25 , 

o., 

860 , 2085 , 

25 , 

o., 


INITIAL GUESS FOR V-VELOCITY (NREC) 
0 , 

DBC FOR V 
35 , 


1 , 

1 . 

1 . 

0 . 

2 , 

2 , 

1 . 

0 . 

3 , 

3 , 

1 , 

0 . 

4 , 

4 , 

1 . 

0 . 

5 , 

5 , 

1 , 

0 . 

6 , 

6 , 

1 , 

0 . 

7 , 

7 , 

1 , 

0 . 

8 , 

8 , 

1 . 

0 . 

9 , 

9 , 

1 . 

0 . 

10 , 

10 , 

1 , 

0 . 

11 . 

11 , 

1 , 

0 . 

12 , 

12 , 

1 , 

0 . 

13 , 

13 , 

1 . 

0 . 

14 , 

14 , 

1 , 

0 . 

15 , 

15 , 

1 , 

0 . 

16 , 

16 , 

1 , 

0 . 

17 , 

17 , 

1 . 

0 . 

18 , 

18 , 

1 . 

0 . 

19 , 

19 , 

1 . 

0 . 

20 , 

20 , 

1 , 

0 . 

21 , 

21 , 

1 , 

0 . 

22 , 

22 , 

1 . 

0 . 

23 , 

23 , 

1 , 

0 . 

24 , 

24 , 

1 , 

0 . 

25 , 

25 , 

1 , 

0 . 

1 . 

501 , 

25 , 

0 . 

25 , 

525 , 

25 , 

0 . 

526 , 

571 , 

15 , 

0 . 

571 , 

585 , 

1 . 

0 . 

586 , 

751 , 

15 , 

0 . 

600 , 

765 , 

15 , 

0 . 

766 , 

779 , 

1 , 

0 . 

780 , 

822 , 

14 , 

0 . 

836 , 2061 , 

25 , 

0 . 

860 , 2085 , 

25 , 

0 . 
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IA04 (PBCDAT, IPN0DE(l-2)) --- 

0., -571, 571, 

END OF INPUT DATA 

****PROCESSOR FOR NAVIER- STOKES EQUATIONS **** 

****END OF RUN **** 

******************************** BOTTOM OF DATA ******************************** 
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APPENDIX III 


DESCRIPTION OF THE SUBROUTINES 

INITAL - Initialize the dimensioned variables. 

BLKDAT - Define the program control parameters, and set the Gauss numerical 
quadrature data in each coordinate direction. 

o DATLIB - Define the flow element to be used, and set the Gauss Numerical 

quadrature data for the computational element. 

PREP - Prepare the input data. 

RNODE - Generate the node coordinate data. 

RELEM - Generate the node connectivity data. 

RINIT - Read in the initial guess. 

RBC1 - Read in the boundary condition data. 

FEMDAT - Read in the re-start data. 

ISOPEL - Compute the interpolation polynomials and the derivatives. 

LSHP1 - Shape functions for one -dimensional linear element. 

LSHP2 - Shape functions for one -dimensional quadratic element. 

SHAP01 - Shape function for two-dimensional constant element. 

SHAP02 - Shape functions for triangular element. 

SHAP23 - Shape functions for bi-quadratic quadrilateral element. 

PROCES - Processor for Navier- Stokes equations. 

PFRONT - Pre-processor for the frontal solver. 

SHPLIB - Save the shape functions on a disk file (logical unit = 2), and 
read the data whenever necessary. 

S1FL0W - Create the sequential degree -of- freedom number for each flow 
variable, and compute the total degrees of freedom. 

■ SEQVFL - Include boundary conditions into the global solution vector. 

SFLOW - Solve the Navier- Stokes equations iteratively. 

FRONTS - Frontal solver. 

ELEMFL - Compute the element system of equations. 
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SCNVFL - Check the convergence . 

SPRS4 - Compute the nodal pressure. 

PFLOW - Print out the computational results. 


USER 


- Load the coordinate data and the flow variables for each element. 
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